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1  Introduction 


Breast  cancer  is  the  most  frequently  diagnosed  malignancy  among  women  in  the  United 
States  [1].  In  1996  the  American  Cancer  Society  estimates  that  184,300  women  will  be 
newly  diagnosed  with  breast  cancer  and  that  44,300  will  die  from  the  disease  [1].  Breast 
cancer  accounts  for  31%  of  all  cancers  detected  and  17%  of  all  cancer  deaths,  and  ranks  as 
the  second  leading  cause  of  death  from  cancer  among  women  in  the  United  States  [1].  Five 
year  survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being 
localized,  falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [2].  The  early 
detection  of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce 
breast  cancer  mortality. 

Mammography’s  role  is  the  early  detection  of  breast  cancer.  Although  more  accurate  than 
any  other  modality,  existing  techniques  for  mammography  only  find  80  to  90%  of  the  breast 
cancers.  Moreover,  in  7  to  10%  of  cases,  the  cancer  will  not  be  visible  on  the  mammogram. 
It  has  been  suggested  that  mammograms  as  normally  viewed,  display  only  about  3%  of  the 
total  information  detected.  Perception  is  a  problem  particularly  for  patients  with  dense 
fibroglandular  patterns.  The  importance  of  diagnosis  of  breast  cancer  at  an  early  stage  is 
critical  to  patient  survival.  The  general  inability  to  detect  small  tumors  and  other  salient 
features  within  mammograms  motivates  our  investigation. 

The  goal  of  this  project  is  to  develop  an  interactive  diagnostic  tool  for  radiologists  that 
shall  refine  the  perception  of  mammographic  features  (including  lesions,  masses  and 
calcifications)  and  improve  the  accuracy  of  diagnosis.  By  improving  the  visualization  of 
breast  pathology  we  can  increase  the  chances  of  early  detection  of  breast  cancers  (improve 
quality)  while  requiring  less  time  to  evaluate  mammograms  for  most  patients  (lower  costs). 
We  are  investigating  a  methodology  for  accomplishing  mammographic  feature  analysis 
through  multiresolution  representations.  We  have  shown  that  overcomplete  representations 
may  be  identified  and  used  to  detect  a,nd  enhance  specific  mammographic  features  within  a 
continuum  of  scale  space.  Such  “focused”  reconstructions  may  complement  existing 
modalities  and  allow  a  radiologist  to  examine  interactively  diagnostic  features  within  a 
selected  scale  space.  Similar  to  traditional  coarse  to  fine  matching  strategies,  the 
radiologist  may  first  choose  to  look  for  coarse  features  (e.g.  dominant  mass)  within  low 
frequency  levels  of  a  wavelet  transform  and  later  examine  finer  features  (e.g. 
microcalcifications)  at  higher  frequency  levels. 

A  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  subtle  difference  in 
x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [3].  This  fact 
makes  the  detection  of  small  malignancies  problematical,  especially  in  younger  women  who 
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have  denser  breast  tissue.  Although  calcifications  have  high  inherent  attenuation 
properties,  their  small  size  also  results  in  a  low  subject  contrast  [4].  As  a  result,  the 
visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be  a  problem 
in  mammography  as  it  is  currently  performed  using  analog  film. 

In  this  report,  we  describe  in  detail  results  accomplished  during  the  past  year.  Presented 
below  in  executive  summary,  we  encapsulate  progress  related  to  the  specific  tasks  and 
objectives  stated  in  Phases  III,  IV  and  V  of  our  Statement  of  Work. 

In  the  sections  below,  we  describe  in  overview,  our  processing  algorithms,  experimental 
methods  and  sample  results  obtained.  In  addition,  we  list  in  summary,  publications  of  our 
researchers  accomplished  during  the  past  year  of  our  investigation.  Finally,  we  include  our 
response  to  the  reviewers  comments  and  suggestions  regarding  contractual  and  technical 
issues  in  reference  to  last  year’s  report. 

1.1  Overview  of  Contents 

Enhancement  via  Fusion  of  Mammographic  Features 

Existing  methods  of  mammographic  image  enhancement  can  be  divided  roughly  into  two 
categories:  (1)  methods  aimed  at  better  visualization  of  all  features  present  in  an  image, 
and  (2)  methods  that  target  specific  features  of  importance  (e.g.,  microcalcifications, 
stellate  lesions).  Methods  from  the  first  category  are  not  optimized  for  a  specific  type  of 
cancer,  while  methods  in  the  second  category  are  typically  tailored  for  a  single  type  of 
malignancy  only.  We  present  a  synthesis  of  the  two  paradigms  by  means  of  image  fusion  to 
overcome  these  shortcomings  and  practical  limitations. 

Both  processing  for  enhancement  of  selected  features  and  fusion  of  the  resultant  images 
were  accomplished  within  a  single  wavelet  transform  framework.  We  chose  a  multiscale 
spline  derivative-based  transform  with  wavelets  equal  to  the  second  derivative  of  a  central 
B-spline  function.  Such  a  transform  not  only  includes  the  two-dimensional  discrete  dyadic 
wavelet  transform  decomposition  as  its  subset,  but  also  enables  incorporation  of  a  variety 
of  known  wavelet  based  image  enhancement  methods.  Enhancement  and  fusion  of 
mammographic  features  within  a  single  transform  framework  enabled  both  efficient 
implementation  and  independent  optimization  of  processing  modules. 

Detection  of  Subtle  Masses  in  Mammography  via  Precision  Multiscale  Analysis 

We  introduce  a  continuous  scale  wavelet  detector  for  identifying  masses  (possible  breast 
cancers)  in  mammograms.  Continuous-scale  wavelet  algorithms  have  been  discussed  in  the 
past,  however  this  is  the  first  reported  algorithm  that  uses  a  scaled  version  of  the  same 
mother  wavelet  at  each  scale  of  analysis.  This  single  mother  wavelet  property  leads  to  a 
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simpler  implementation  and  a  more  direct  application  of  detection  theory  to  recognition 
problems  than  traditional  multiscale  analysis.  In  addition,  we  show  that  a  continuous-scale 
search  is  necessary  for  computer  aided  diagnosis  of  mammography  since  traditional 
solutions  using  dyadic  scales  (powers  of  two)  either  fail  to  detect  some  masses  or  signal  too 
many  false  alarms.  Our  novel  wavelet  detector  combines  a  wavelet  formulation  with  the 
classical  theory  of  constant  false  alarm  rate  (CFAR)  detectors.  We  show  that  our  algorithm 
is  able  to  detect  masses  in  dense  mammograms  that  could  not  be  seen  using  conventional 
windowing  and  leveling  or  other  traditional  methods  of  contrast  enhancement. 

Coherence  of  Multiscale  Features  for  Digital  Mammogram  Enhancement 

Mammographic  features  are  extracted  according  to  their  coherence  and  orientation.  In  this 
chapter  of  the  report,  an  artifact  free  enhancement  algorithm  based  on  overcomplete 
multiscale  wavelet  analysis  is  presented.  First,  an  image  was  decomposed  using  a  fast 
wavelet  transform  algorithm.  At  the  same  time,  the  energy  and  phase  information  at  each 
level  are  determined  using  a  set  of  separable  steerable  filters.  Then,  a  measure  of  coherence 
within  each  level  is  obtained  by  weighting  an  energy  measure  with  the  ratio  of  projections 
of  the  energy  within  a  specified  window  onto  the  central  point  of  a  window  with  respect  to 
the  total  energy  within  each  window.  Finally,  a  nonlinear  operation,  integrating  coherence 
and  orientation  information,  is  applied  to  modify  transform  coefficients  within  distinct 
levels  of  analysis.  These  modified  coefficients  are  then  reconstructed,  via  an  inverse  fast 
wavelet  transform,  resulting  in  an  improved  visualization  of  mammographic  features.  The 
novelty  of  this  algorithm  lies  in  its  detection  of  directional  features  and  reconstruction 
without  artifacts.  This  algorithm  represents  the  “state-of-the-art”  of  multiscale 
enhancement  approaches,  and  is  the  algorithm  of  choice  for  our  ongoing  ROC  studies. 

A  Lifting  Algorithm  for  Overcomplete  Wavelet  Representations  for  Interactive 
Feature  Analysis 

The  Fourier  transform  has  been  the  traditional  mathematical  tool  in  wavelet  processing 
because  operations  of  a  wavelet  family  defined  by  translates  and  dilates  become  algebraic 
operations  in  Fourier  space.  However,  the  Fourier  transform  is  not  efficient  if  the  filter 
length  involved  is  fairly  short,  and  may  not  be  applicable  to  problems  on  more  complex 
domains.  In  addition,  since  wavelet  functions  dealing  with  bounded  inputs  are  not 
translation  invariant,  the  Fourier  transform  is  best  not  used  as  a  construction  tool  of 
wavelets.  To  overcome  these  drawbacks,  we  propose  a  flexible  and  efficient  tool  of  analysis, 
namely  overcomplete  lifting  scheme,  which  will  be  used  to  compute  overcomplete  wavelet 
representations  (in  real  time)  for  interactive  mammographic  feature  analysis. 
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1.2  Publications 

Below,  we  list  in  summary,  an  update  of  publications  accomplished  during  1997. 

1.  D.  Chen,  A.  Laine,  J.  G.  Harris,  and  W.  Huda,  “A  continuous  scale  discrete  wavelet 
transform  for  the  detection  of  spicular  masses  in  mammograms,”  submitted  to  IEEE 
Transactions  on  Signal  Processing,  February  1997. 

2.  D.  loannou,  W.  Huda,  A.  F.  Laine,  and  I.  Koren,  “Enhancement  of  computer  simulated 
masses  for  mammography  via  wavelet  analysis,”  submitted  to  IEEE  Transactions  on 
Medical  Imaging,  January  1997. 

3.  I.  Koren,  A.  Laine,  and  F.  Taylor,  “Enhancement  via  fusion  of  mammographic  features, 
submitted  to  IEEE  International  Conference  on  Image  Processing,  Chicago,  IL,  October 
1998. 

4.  M.  Shim  and  A.  Laine,  “Overcomplete  lifted  wavelet  representations  for  multiscale 
feature  analysis,”  submitted  to  IEEE  International  Conference  on  Image  Processing, 
Chicago,  IL,  October  1998. 

5.  1.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  M.  Akay  (Editor),  Time-Prequency  and  Wavelets  in  Biomedical  Signal 
Engineering^  New  York,  NY:  IEEE  Press,  1998,  pp.  425  449. 

6.  S.  Schuler  and  A.  Laine,  “Hexagonal  QMF  banks  and  wavelets,”  in  M.  Akay  (Editor), 
Time-Frequency  and  Wavelets  in  Biomedical  Signal  Engineering,  New  York,  NY:  IEEE 
Press,  1998,  451-472. 

7.  D.  Chen,  C.-M.  Chang,  and  A.  Laine,  “Detection  and  enhancement  of  small  masses  via 
precision  multiscale  analysis,”  in  Proceedings  of  the  Third  Asian  Conference  on  Computer 
Vision,  Hong  Kong,  PRC,  January  1998,  vol.  1,  pp.  192-199. 

8.  A.  Laine,  W.  Huda,  J.  Honeyman,  and  B.  Steinbach,  “Wavelet  representations  for  digital 
mammography,”  in  Proceedings  of  the  Department  of  Defense  Breast  Cancer  Research 
Program  Meeting  Era  of  Hope,  Washington,  D.C.,  October-November  1997,  vol.  1,  pp. 
105-106. 

9.  1.  Koren,  A.  Laine,  and  F.  Taylor,  “An  overcomplete  enhancement  of  digital 
mammograms,”  in  Proceedings  of  the  Department  of  Defense  Breast  Cancer  Research 
Program  Meeting  Era  of  Hope,  Washington,  D.C.,  October-November  1997,  vol.  1,  pp. 
107-108. 

10.  X.  Zong,  A.  Meyer-Baese,  and  A.  Laine,  “Multiscale  segmentation  through  a  radial 
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basis  neural  network,”  in  Proceedings  of  the  IEEE  International  Conference  on  Image 
Processing,  Santa  Barbara,  October  1997,  vol.  3,  pp.  400-403. 

11.  C.-M.  Chang  and  A.  Laine,  “Enhancement  of  mammograms  from  oriented 
information,”  in  Proceedings  of  the  IEEE  International  Conference  on  Image  Processing, 
Santa  Barbara,  October  1997,  vol.  3,  pp.  524-527. 

1.3  Responses  to  Technical  and  Contractual  Issues 

With  regard  to  the  ROC  studies  mentioned  in  Phase  V  of  our  SOW,  we  include  below  a 
description  of  an  ongoing  investigation.  The  ROC  studies  have  been  delayed  due  to  the 
fact  that  the  principal  investigator  and  two  of  the  three  principal  collaborators  of  the 
project  within  the  Department  of  Radiology  at  the  University  of  Florida,  resigned  their 
positions  to  accept  new  academic  positions  at  other  institutions.  The  length  of  the  project, 
has  been  extended  to  accommodate  this  adjustment  and  allow  successful  completion  of  the 
project.  In  addition,  new  personnel  affiliated  with  Columbia-Presbyterian  Hospital, 
Department  of  Radiology,  Columbia  University,  will  assist  in  carrying  out  the  final  ROC 
studies  described  in  Phase  V  of  the  Statement  of  Work.  We  expect  that  these  studies  will 
be  complete  within  six  months  time. 

Performance  Study:  We  will  use  an  existing  national  mammography  database 
(currently  under  development  at  the  University  of  South  Florida  and  Massachusetts 
General  Hospital  under  the  supervision  of  Dr.  Bowyer  and  Dr.  Kopans,  respectively)  in 
order  to  carry  out  a  receiver  operating  characteristic  (ROC)  study  to  evaluate  the  effect  of 
multiscale  enhancement  (MSE)  on  the  diagnostic  performance  of  radiologists.  The  study 
will  use  360  four- view  clinical  cases:  180  cases  with  cancer  proved  by  means  of  biopsy  and 
180  cases  with  negative  findings  at  examination  and  follow-up.  To  measure  the  radiologists 
performance  with  and  without  MSE,  the  360  cases  will  be  divided  into  three  groups  of  120 
cases  each  (randomly  chosen.)  The  study  will  be  conducted  in  two  phases.  In  each  phase 
all  360  cases  will  be  prepared  on  an  alternator  in  random  order.  During  the  first  phase, 
cases  in  the  first  group  will  be  presented  with  four  original  views  and  enhanced  images, 
while  cases  in  the  second  and  third  groups  will  be  presented  with  only  four  original  views. 
Similarly,  during  the  second  phase,  cases  in  the  second  group  will  be  presented  with  four 
original  views  and  MSE  images,  while  cases  in  the  first  and  third  groups  will  be  presented 
with  only  four  original  views.  Eight  weeks  will  be  allowed  between  the  two  phases  of  the 
study  to  minimize  any  potential  learning  by  the  observers.  Each  phase  will  be  divided  into 
six  reading  sessions  on  separate  days.  During  each  session  an  observer  will  read  one-sixth 
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Table  1:  Case  distribution. 


Total  number  of  cases 

360 

Number  of  normal  cases 

180 

Number  of  malignant  cases 

180 

Number  of  phases 

2 

Number  of  groups  per  phase 

3 

Number  of  cases  per  group 

120 

Number  of  sessions  per  phase 

6 

Number  of  cases  per  session 

60 

Table  2:  Reading  conditions  for  each  group  of  cases. 


Group 

Cases  with 

MSE 

Phase  1 

Phase  2 

1 

All 

None 

2 

None 

All 

3 

None 

None 

of  the  mammograms  (60  cases)  sequentially.  Table  1  and  Table  2  show  a  summarized 
description  of  the  performance  study. 

Six  mammographers  will  participate  in  the  ROC  study.  None  of  the  observers  will  be 
informed  of  the  prevalence  rate  or  have  previous  exposure  to  the  cases  in  the  study.  In  each 
session  the  observers  will  be  required  to  rate  for  each  case  their  confidence  regarding  the 
presence  of  malignancy  using  a  five  category  scale  (1  =  definitely  not  present,  2  =  probably 
not  present,  3  =  possibly  present,  4  =  probably  present,  5  =  definitely  present.) 

ROC  analysis  [5,  6]  will  be  employed  to  compare  the  observers  performance  in  detecting 
the  presence  of  malignancy  with  and  without  MSE.  For  each  observer,  the  accuracy  index 
will  be  calculated  for  each  group  of  cases.  The  differences  in  A,  between  phases  one  and 
two  for  a  given  group  will  indicate  if  the  performance  of  an  observer  increased  or  decreased 
with  MSE.  Table  2  summarizes  the  reading  conditions  for  each  group  of  cases.  The  first 
group  will  measure  the  change  in  A^  when  the  cases  are  first  presented  with  the 
computerized  outputs  as  ancillary  films.  The  second  group  will  measure  the  change  in  Az 
when  the  cases  are  first  presented  without  the  computerized  outputs.  The  third  group 
serves  as  a  control  and  will  measure  the  change  (if  any)  when  the  cases  are  presented  to  the 
observers  twice  without  computerized  outputs.  By  means  of  this  study  we  will  learn 
whether  or  not  mammographers  perform  better  in  the  detection  of  small  malignancies 
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when  MSE  images  are  made  available. 


2 


Body 

2.1  Enhancement  via  Fusion  of  Mammographic  Features 

2.1.1  Introduction 

In  mammography,  early  detection  of  breast  cancer  relies  upon  the  ability  to  distinguish 
between  malignant  and  benign  mammographic  features.  The  detection  of  small 
malignancies  and  subtle  lesions  is  often  difficult.  Contrast  enhancement  can  make  more 
obvious  unseen  or  barely  seen  features  of  a  mammogram  without  requiring  additional 
radiation. 

Existing  methods  of  mammographic  image  enhancement  can  be  divided  roughly  into  two 
categories:  (1)  methods  aimed  at  better  visualization  of  all  features  present  in  an  image 
[7,  8,  9,  10],  and  (2)  methods  that  target  specific  features  of  importance  (e.g., 
microcalcifications  [11,  12,  13],  stellate  lesions  [14]).  Methods  from  the  first  category  are 
not  optimized  for  a  specific  type  of  cancer,  while  methods  in  the  second  category  are 
typically  tailored  for  a  single  type  of  malignancy  only.  In  this  chapter,  we  present  a  novel 
approach  which  overcomes  these  shortcomings  and  problematic  limitations. 

Since  diagnostic  features  in  mammograms  appear  in  a  variety  of  shapes  and  sizes, 
traditional  image  enhancement  techniques  such  as  histogram  equalization  and  unsharp 
masking  seldom  produce  satisfactory  results  and  are  clearly  outperformed  by  more 
sophisticated  methods  [9,  10].  Due  to  their  ability  to  analyze  images  across  different  scales, 
wavelet-based  techniques  have  become  particularly  popular  for  processing  of  mammograms 
[7,  8,  9,  11,  12,  13].  The  majority  of  researchers  lean  toward  redundant  wavelet 
representations  which  do  not  introduce  significant  processing  artifacts  such  as  aliasing 
inherent  to  orthogonal  and  biorthogonal  wavelet  transforms.  Here,  we  continue  with  a 
generalization  of  the  discrete  dyadic  wavelet  transform  started  in  our  previous  report  and 
concentrate  on  wavelets  that  are  equal  to  the  second  derivative  of  a  central  B-spline 
function. 

This  chapter  is  organized  as  follows:  in  Sections  2.1.2  and  2.1.3  a  discrete  dyadic  wavelet 
transform  in  one  and  two  dimensions,  respectively,  is  presented.  Next,  Section  2.1.4 
explains  a  novel  approach  of  transform-based  feature  enhancement  via  image  fusion. 

Section  2.1.5  then  describes  transform  implementation  in  detail.  Finally,  Section  2.1.6 
provides  a  conclusion. 

2.1.2  1-D  Discrete  Dyadic  Wavelet  Transform  Revisited 

Let  us  begin  with  a  brief  review  of  properties  of  the  one-dimensional  discrete  dyadic 
wavelet  transform  as  described  in  our  previous  report,  but  included  here  for  completeness. 
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The  dyadic  wavelet  transform  of  a  function  s(a;)  G  L‘^{R)  is  defined  as  a  sequence  of 
functions 


(1) 

where  Wms{x)  =  s*  '0m(a:)  =  s{t)  'ipmix-t)  dt,  ^ra{x)  =  is  a  wavelet 

expanded  by  a  dilation  parameter  (or  scale)  2"”,  and  denotes  the  convolution  operator. 
To  ensure  coverage  of  the  frequency  axis  the  requirement  on  the  Fourier  transform  of 
ipm{x)  is  the  existence  of  >  0  and  <  oo  such  that 

OO 

m=-oo 

is  satisfied  almost  everywhere.  The  constraint  on  the  Fourier  transform  of  the  (nonunique) 
reconstructing  function  x(x)  is 

OO 

«2’»a,)  =  1, 

m=—oo 

A  function  s{x)  can  then  be  completely  reconstructed  from  its  dyadic  wavelet  transform 
using  the  identity 

OO 

s{x)  =  ^  ^  bFmS  *  Xm{^\ 

m=—oo 

where  Xm{x)  =  2~'^x{^~'^^)- 

In  numerical  applications,  processing  is  performed  on  discrete  rather  than  continuous 
functions.  When  the  function  to  be  transformed  is  in  the  discrete  form,  the  scale  2”^  can  no 
longer  vary  over  all  m  G  Z.  Finite  sampling  rate  prohibits  the  scale  from  being  arbitrarily 
small,  while  computational  resources  restrict  the  use  of  an  arbitrarily  large  scale.  Let  the 
finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  2^. 

The  smoothing  of  a  function  s{x)  G  L‘^{R)  is  defined  as 

Sm^i^X^  S  *  (f)ffi{x^  ^ 

where  ^m{x)  =  2~'^(f){2~'^x)  with  me  Z,  and  (j){x)  is  a  smoothing  function  (i.e.,  its 
integral  is  equal  to  1  and  (j){x)  0  as  |a;|  ^  oo). 

Mallat  and  Zhong  [15]  selected  a  real  smoothing  function  (/){x),  whose  Fonrier  transform 
satisfied 

OO 

|,AMP  =  5;^(2’"a.)x(2’"a;).  (2) 

m—1 
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In  addition,  it  was  shown  that  any  discrete  function  of  finite  energy  s(n)  G  can  be 

written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1,  i.e.,  s{n)  =  Sof{n), 
where  f{x)  G  L‘^{R)  is  not  unique.  Thus,  the  discrete  dyadic  wavelet  transform  of  s{n)  for 
any  coarse  scale  2^  was  defined  as  a  sequence  of  discrete  functions 

{SMf{n  +  s),  {Wmf{n  +  s)}me[x,M]}neZ^ 


where  s  is  a  ^^{x)  dependent  sampling  shift. 

The  above  initialization  s(ri)  =  Sof{n)  is  rather  standard  in  the  discrete  wavelet  transform 
computation  [16],  although  it  yields  correct  results  (i.e.,  the  discrete  wavelet  transform  is 
equal  to  the  samples  of  its  continuous  counterpart)  only  when  s{n)  =  Sos{n).  Here,  we  will 
concentrate  on  wavelets  which  are  derivatives  of  spline  functions  and  this  will  lead  us  to  a 
simple  initialization  procedure  [17]  that  alleviates  the  above  problem. 

For  a  certain  choice  of  wavelets,  the  discrete  dyadic  wavelet  transform  can  be  implemented 
within  a  fast  hierarchical  digital  filtering  scheme.  Next,  we  shall  summarize  the  relations 
between  filters,  wavelets,  and  smoothing  functions. 

First,  let  us  introduce  a  real  smoothing  function  ip{x)  such  that  (2)  can  be  rewritten  as 

OO 

^{cv)  ip{u})  =  ^  '^{2'^u)  x(2™'a;),  (3) 

m—0 


and  let  us  set 


(j>{x)  =  Pp{x)  =  ^  + 


P+1 


—  i  \  u\ X  + 


p+1 


—  i 


(4) 


where  Ppix)  denotes  a  central  B-spline  of  order  p.  With  the  choice  (4),  we  restrict  ourselves 
to  wavelets  which  are  spline  functions. 

Computing  (3)  for  the  finest  two  scales  shows  that 


-iplu)  x{lo)  =  Pp{uj)  (p{uj)  -  Pp{2oj)  (p{2u). 

Pp{2uj)  can  be  related  to  ^p{u))  by  expressing  I3p{2uj)  as  (cf.  Proposition  1  of  [17]) 


(5) 


$p{2co)  = 


sin(a;) 

^  Isin  (f) 


p  +  1 


Sin 


(!) 


P  +  1 


2 


and  using  E“.. 


jdp(2u)  =  (cos  (^))^  ++)■ 


(6) 


^Note  that  the  sum  index  determines  the  range  of  scales  of  the  discrete  transform:  using  (2)  we  have 
ifi(2LL>)  and  x(2u;)  at  the  finest  scale  of  the  transform,  while  for  (3)  we  get  and  x(w). 
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(Note  that  a  relation  similar  to  (6)  can  be  derived  for  integer  scales  provided  that  the 


dilation  parameter  and  the  order  p  are  not  both  even  [17].) 

Let  F(u)  be  a  digital  filter  frequency  response  and  let  Fs(u>)  =  e^^^F(Lo). 

If  we  choose 

il)(u})  =  G-s(<^)  ^p(<x), 

(7) 

(p(2u))  =  Ls(ijo)(p{uj), 

(8) 

X(ix)  =  Ks(lo)(p(uj), 

(9) 

and 

BM  =  e>"(cos(|))’’’‘‘. 

(10) 

where  s  G  {0,  is  a  filter  dependent  sampling  shift  needed  for  g{n),  l{n),  k{n),  and  h{n) 
to  be  FIR  filters,  and  insert  Equations  (6)-(10)  into  (5),  we  observe  the  relation  between 
frequency  responses  of  the  filters 


G{co)K{u;)  +  H{u)L{u:)  =  1.  (11) 

Similar  to  orthogonal  and  biorthogonal  discrete  wavelet  transforms,  the  discrete  dyadic 
wavelet  transform  can  be  implemented  within  a  hierarchical  filtering  scheme.  To  derive 
such  a  digital  filtering  scheme,  let  us  assume  that  s{uj)  from  (1)  is  bandlimited  to  [-tt,  tt]. 
Using  Shannon’s  sampling  theorem  [18]  and  (7)  in  the  definition  of  the  dyadic  wavelet 
transform  (1)  with  m  =  0  shows 

/OO  ^  ^ 

^  s(i)sinc(t  -  i)  ^  5'_5(m)/?p(a;  -  t  -  m)  dt 

OO  m=— OO 

By  making  use  of  the  fact  that  the  cardinal  spline  functions  rjr  (x)  tend  to  the  sine  function 
as  their  order  r  approaches  infinity  [19],  and  employing 

OO 

Vr(x)  =  '^  b~^(i)/3r(x  -  i),  (12) 

i=—oo 

where  b~^(n)  denotes  a  direct  B-spline  filter  [20],  we  can  write 

Wqs(u))  ~  S(io)B~^(uj)^r(<^)  ^p(uj)G-s(^), 
or,  by  using  (6)  and  (10), 

m— 1 

^  S(u<)  Bp+.«(a,)  G_.(2’"c.)  J]  (13) 

i=0 
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Figure  1:  Filter  bank  implementation  of  a  one-dimensional  discrete  dyadic  wavelet  transform 
decomposition  (left)  and  reconstruction  (right)  for  three  levels  of  analysis. 


Equation  (13)  entirely  specifies  the  discrete  dyadic  wavelet  transform  decomposition,  while 
the  reconstruction  follows  from  (5)-(10).  Three  levels  of  a  filter  bank  implementation  are 
shown  in  Figure  1.  Note  that  the  initialization  is  the  same  as  the  one  proposed  in  [17]  and 
that  except  for  the  prefiltering  and  postfiltering,  this  scheme  is  implementing  “algorithme  a 
trous”  [21].  Noninteger  shifts  at  scale  1  for  filters  with  s  =  |  are  rounded  to  the  nearest 
integer. 

Next,  we  choose  filters  in  the  filter  bank  implementation  of  the  discrete  dyadic  wavelet 
transform.  As  previously  mentioned,  we  are  interested  in  wavelets  which  are  derivatives  of 
spline  functions.  From  the  property  of  the  central  B-spline  functions  [22] 

d(5p{x)  ( r.  —  R  .  /^.r-  — 


it  follows  that  G{ijo)  in  (7)  is  the  Fourier  transform  of  the  difference  operator  centered 
around  —s: 


G{ijj)=  ( 2j  sin 


where  d  is  the  order  of  the  derivative  and  the  sampling  shift  for  this  filter  is  s  —  ^  ■ 

Since  H{ijo)  was  already  given  by  (10),  the  remaining  two  filters  to  be  determined  are  L{uj) 
and  K{u)).  Both  of  them  are  (as  is  true  for  ip{x)  and  x(a:))  nonunique. 

If  we  choose  L{ijo)  such  that  we  can  express  K{u})  in  terms  of  a  finite  geometric  series 
having  the  smallest  number  of  elements  for  an  arbitrary  p,  we  get 

I  I 

/  1^1  \  /  /a;\\  , 

i(^) = y:  i-ir*'  M  M  (cos  ( x) )  (15 


(p+l){2m-l) 


K(uj) 


e  sin 


\  \  (imod2 


t(-(l))^ 
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Table  3:  Impulse  responses  h{n),  g{n),  l{n),  and  A:(n)  for  p  =  0  and  d  g  {1,2,3}. 


1 

d  =  2 

d  =  3 

n 

h{n) 

l{n) 

k{n) 

di'IT') 

l{n) 

k{n) 

9{n) 

l{n) 

k{n) 

-2 

1 

-1 

0.5 

1 

1 

-3 

-0.125 

0 

0.5 

-1 

0.5 

-0.25 

-2 

0.5 

-0.25 

3 

0.625 

0.0625 

1 

0.5 

0.25 

1 

0.5 

-1 

0.625 

-0.0625 

2 

-0.125 

where  [x\  denotes  the  largest  integer  smaller  than  x,  the  sampling  shift  for  L{u)  is  the 
same  as  the  one  for  H{uj)  (i.e.,  and  the  sampling  shift  for  K{uj)  is  the  same 

as  the  one  for  G{oj). 

Note  that  Equations  (10)  and  (14)-(16)  work  fine  from  the  mathematical  point  of  view, 
but  in  practice  the  reconstruction  may  become  cumbersome  when  both  p  and  d  are  large. 
(The  lengths  of  impulse  responses  h{n),  g{n),  l{n),  and  k{n)  are  p+2,  d+1, 
(p+l)(d-(d+l)mod2)+l,  and  pd+(p+l)(dmod2)+l,  respectively,  while  for  the  frequency 
responses  of  the  decomposition  filters  we  observe  that  limp^oo  \H-s{>^)  \  =  dy,{uj  +  2n7r)  and 
limd^oo(2i)“^  |G'-5(‘^)I  =  ^^i(^  +  (2n+l)7r)  with  neZ.) 

It  is  also  worth  noting  that  K{u)  is  a  lowpass  filter  when  p  is  even  (i.e.,  the  reconstruction 
function  xi^)  is  a  wavelet  only  for  p  odd). 

Tables  3,  4,  and  5  list  impulse  responses  of  the  four  filters  for  p  e  {0, 1, 2}  and  d  e  {1, 2, 3}, 
while  Figure  2  shows  the  corresponding  wavelets  il;{x)  =  Wavelets  from  this 

family  have  a  support  of  length  d+p+1,  regularity  order  p,  and  are  either  symmetric  or 
antisymmetric. 

2.1.3  2-D  Discrete  Dyadic  Wavelet  Transform  Revisited 

In  two  dimensions,  the  dyadic  wavelet  transform  of  a  function  s(x,y)  G  is  defined 

as  a  set  of  functions  [15] 

{W^s{x,y),W^s{x,y)}^^2:^  (17) 

where  Wls{x,  y)  =  s*  ^^(x,  y)  for  i  -  {1,  2}  and  ^/;^(x,  y)  =  2-2"^t/)*(2-”»x,  2-^y)  are 
wavelets  #(x,  ?/)  expanded  by  a  dilation  parameter  2^. 

To  ensure  coverage  of  the  frequency  space  there  must  exist  A2  >  0  and  B2  <  00  such  that 

00  2 

^2  <  E  ^2 
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Table  4;  Impulse  responses  h 


n 


),  g{n),  l{n),  and  k{n)  for  p=l  and  d  G  {1, 2, 3} 


n 

h{n) 

l{n) 

d  =  l 

d  =  2 

9{n) 

k{n) 

9in) 

k{n) 

-1 

0 

1 

2 

0.25 

0.5 

0.25 

1 

-1 

-0.0625 

-0.3125 

0.3125 

0.0625 

1 

-2 

1 

-0.0625 

-0.375 

-0.0625 

n 

d  =  3 

h{n) 

9{1T') 

l{n) 

k{n) 

-3 

-2 

-1 

0 

1 

2 

3 

0.25 

0.5 

0.25 

1 

-3 

3 

-1 

-0.015625 

-0.09375 

0.265625 

0.6875 

0.265625 

-0.09375 

-0.015625 

0.00390625 

0.04296875 

0.1015625 

-0.1015625 

-0.04296875 

-0.00390625 

Table  5: 


.  XiJ.1 

d  = 

1 

d  = 

2 

n 

h{n) 

9{n) 

l{n) 

k{n) 

9{n) 

l{n) 

k{n) 

-2 

-1 

0.125 

0.375 

1 

0.125 

-0.015625 

-0.109375 

1 

0.125 

-0.015625 

-0.125 

0 

0.375 

-1 

0.375 

-0.34375 

-2 

0.375 

-0.46875 

1 

0.125 

0.375 

0.34375 

1 

0.375 

-0.125 

2 

3 

0.125 

0.109375 

0.015625 

0.125 

-0.015625 

d  =  3 

n 

h{n) 

9{n) 

l{n) 

k{n) 

-4 

-0.001953125 

0.000244140625 

-3 

-0.017578125 

0.003662109375 

-2 

0.125 

1 

-0.0703125 

0.0263671875 

-1 

0.375 

-3 

0.085937 

0.0908203125 

0 

0.375 

3 

0.50390625 

0.13037109375 

1 

0.125 

-1 

0.50390625 

-0.13037109375 

2 

0.0859375 

-0.0908203125 

3 

-0.0703125 

-0.0263671875 

4 

-0.017578125 

-0.003662109375 

5 

-0.001953125 

-0.000244140625 
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k 


is  satisfied  almost  everywhere.  If  (nonunique)  functions  x\^,y)  and  x^{x,y)  am  chosen 
such  that  their  Fourier  transforms  satisfy 


oo  2 

m=— oo 

the  function  s{x,y)  may  be  reconstructed  from  its  dyadic  wavelet  transform  by 

OO  2 

s{x,y)=  ^  *  Xm > 


where  Xmi^^y)  =  2  "^y). 

However,  when  processing  discrete  functions  the  scale  2-  may  no  longer  vary  over  all 
meZ.  Let  the  finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  be  2  .  Let  us 
similar  to  [15],  introduce  a  real  smoothing  function  (t>{x,y)  such  that  its  Fourier  transform 

satisfies 


OO  2 

m=0  i=l 


Un 


,2™u;,)x'(2"^a;„2-a;,). 


(18) 


Here,  as  in  one  dimension,  a  finite  energy  discrete  function  {s{n,,ny)  e  P(^'))  can  be 
written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1:  s(n,,  Uy) 

=  5o/(n^,%),  where  f{x,y)  G  is  not  unique,  and  5^/(x,2/)  =  /  *  <t>m{x,y).  This 

led  Mallat  and  Zhong  [15]  to  a  two-dimensional  analog  of  the  one-dimensional  definition  of 
the  discrete  dyadic  wavelet  transform: 


{Sm—iJ  iV'x 


+  5 


),  {W^/(nj,  +  s,  nj,  +  s),  Wlf{n^  +  s,  %  -f  s)}^e[o,M-i]}  (^Tlx  jRy  )  e 


We  will  use,  as  in  Section  2.1.2,  a  spline-based  initialization  procedure. 

To  implement  a  multidimensional  discrete  dyadic  wavelet  transform  within  a  fast 
hierarchical  digital  filtering  scheme,  the  wavelets  were  chosen  to  be  separable  products  o 

one-dimensional  functions: 

^\x,y)  =  Hx)  <l>{y)^ 

'tp'^{x,y)  =  'ip{y)  (t>{x), 

where  (/)(a;)  and  ^{x)  were  chosen  as  described  in  Section  2.1.2  (i.e.,  (t>{x)  =  Pp{x)  and 
specified  by  (7)). 

2Xs  in  Section  2.1.2,  we  put  the  finest  scale  of  the  transform  at  m  =  0. 
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Table  6:  Impulse  responses  ti{n)  for  p  G  {0, 1,  2}. 


n 

p=0 

P  =  1 

p=2 

-3 

0.0078125 

-2 

0.03125 

0.046875 

-1 

0.125 

0.125 

0.1171875 

0 

0.75 

0.6875 

0.65625 

1 

0.125 

0.125 

0.1171875 

2 

0.03125 

0.046875 

3 

0.0078125 

From  (19),  (20),  and  (7),  we  may  write 

(21) 

"0  i^^x^^y')  s (^2/) /^p(^a;)/^p (^j/) )  (22) 

where  G(uj)  is  given  by  (14)  for  d  G  {1,2}.  Choosing 

X  i^X^^y)  ~  ^^s{^x)Tl(iL)y')Pp(u>x)Pp{u)y),  (23) 

X  ij^xi^y)  ~  Kg(^(jJy')Ti(^UJx')/3p{LOx}0p{ix)y')i  (24) 

where  K{oj)  and  Ti{u>)  are  digital  filter  frequency  responses,  we  may  compute  (18)  for  the 
finest  two  scales  by 
2 

'^'ip\2Ux,2cOy)  X^{2c0x,‘2^y)  =  |0(wx,a;j/)p  -  \^{2ijJx,2ujy)\^ .  (25) 

i=l 

Inserting  the  terms  defined  by  (21),  (22),  (23),  (24),  (6),  and  (10)  with 
^{uJx,ijJy)  =  0p{oJx)0p{(x)y)  into  (25)  results  in 

K{u}x)G{ux)Ti{u)y)  +  K{LOy)G{ojy)Ti(u)x)  +  \H(ux)\^\H{ujy)\^  =  1.  (26) 


Equation  (26)  represents  the  relation  between  the  frequency  responses  of  the  digital  filters 
used  to  implement  a  multidimensional  discrete  dyadic  wavelet  transform  and  is  a 
multidimensional  analog  to  (11). 

Solving  (26)  for  Ti{u>)  by  substituting  K{lo)G{lo)  from  (11)  yields  the  closed  formula  [15] 

=  +  (27) 

In  Table  6,  we  provide  the  filter  coefficients  for  Ti{u>)  from  (27)  for  p  G  {0, 1, 2}.  All  other 
filters  from  (26)  were  already  specified  in  Section  2.1.2. 

As  in  the  one-dimensional  case,  a  two-dimensional  discrete  dyadic  wavelet  transform  can  be 
implemented  as  a  fast  hierarchical  filtering  scheme.  To  derive  such  an  implementation,  we. 
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similar  to  the  one-dimensional  case  from  Section  2.1.2,  use  the  definition  of  the 
two-dimensional  dyadic  wavelet  transform  (17)  and  require  s(a;^,a;y)  =  0  for  |a;a:l  >  tt  or 
\ujy\  >  TT.  Using  Shannon’s  sampling  theorem  in  two  dimensions  [23],  (21),  (22),  and  m  =  0, 

we  get 

poo  poo 

Wqs{x^  y)=  /  s\nc{tx  —  ix)  sinc(ty  —  %)• 

J  —OO  j  —oo  7*  —  _ry-)  =  — oo 


OO 

g_s{'m)/Sp{x  -  4  -  m)pp{y  -  ty)  dt^  dty 


poo  poo  .  .  /  .  \ 

W^s{x,y)=  /  /  2-/  s{ix,iy)sinc{t3;-ix)sinc{ty-iy)- 

J  —OO  J  —oo 

OO 

•/3p{x  —  tx)  X  9-s{'<xi)ldp{y  —  ty  —  m)  dtx  dty. 

m^—oo 

As  in  one  dimension,  we  make  use  of  the  fact  that  the  cardinal  spline  functions  converge  to 
the  sine  function  as  their  order  r  tends  to  infinity  and  write 

^siUx.l^y)  -  S{uJx,UJy)B^'^{uJx)B;^{uJy)^p+r+l{^x)h+T+\MG-s{.^x) 

and 

\^s{0Jx,U)y)  ^  S{tOx,(^y)B;\0Jx)B;\u}y)^p+r+lMPp+r+lM  G-s{<Xy), 
or  for  me  [0,  M)  and  discrete  signal  processing 

T{W^^s{x,y)Y^  _  }  ^  S{ujx,ujy)  B;\u;x)  B;\u;y)  Bp+r+iM- 


■fit— X 

G-.s{2  tOx)  n  H^,{2^UJx)H-siP‘"^y) 


J^{W^s{x,  y)\x,^n^^y=ny^  "  S{uJx,<^y)  B.^  ^(u;^;)  ^{x>y)  Bp+rJrl{^x)' 

m—1 

■Bp+r+lMG^si^y^^y) 

i=0 

Equations  (28)  and  (29)  describe  the  decomposition  part  of  the  filter  bank  implementation 
of  a  two-dimensional  discrete  dyadic  wavelet  transform.  The  reconstruction  part  can  be 
obtained  from  (23)-(25)  with  (6)  and  (10).  The  entire  filter  bank  implementation  of  the 
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(a) 


(b) 


Figure  3:  Filter  bank  implementation  of  a  two-dimensional  discrete  dyadic  wavelet  transform 
(a)  decomposition  and  (b)  reconstruction  for  two  levels  of  analysis.  denotes  the 

complex  conjugate  of 
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transform  is  shown  in  Figure  3.  Except  for  the  prefiltering  and  postfiltering,  we  readily 
recognize  the  implementation  proposed  in  [15]. 

Using  the  fact  that  a  wavelet  'ipix)  is  equal  to  a  d-th  derivative  of  a  spline  function  Pp+dix), 
(19)  and  (20)  may  be  rewritten  as 


dyd  ? 


where 

d^{x,y)  =  Pp+d{x)  Pp{y) 

and 

d'^ix,y)  =  Pp{x)  Pp+d{y)- 

Let  us  denote  Wjns{x^  y)  =  (W^s{x,  y),  W^s{x,  y)),  V  =  (^,  ^),  A  =  V  —  (^  +  Qy2), 
and  assume  that  ?9*(a:,  y, )  can  be  approximated  by  'd{x,  y)  for  i,  d  G  {1,  2}. 

For  d  =  l  it  then  follows  that 

Wrns{x,  y)  =  2”^V(s  *  dm){x,  y).  (30) 

Thus  for  d=2  we  can  write 

2 

'^W!^s{x,y)  =  2‘^'^/\{s*'dm){x,y).  (31) 

i=l 

— ¥ 

With  ^{x,  y)  being  a  Gaussian,  finding  local  extrema  of  (30)  in  the  direction  of  gradient  V 
corresponds  to  the  filtering  stage  of  a  Canny  edge  detector  [24] ,  and  finding  zero-crossings 
of  (31)  corresponds  to  the  filtering  carried  out  with  a  Marr-Hildreth  edge  detector 
(Laplacian  of  Gaussian)  [25].  (Note  that  both  edge  detectors  involve  postprocessing).  Edge 
detection  based  on  finding  local  extrema  of  W^s(a;,y)  or  zero-crossings  of  W^mS(a;,  y) 
is  therefore  an  approximation  to  the  Canny  or  the  Marr-Hildreth  edge  detector  over  a 
range  of  dyadic  scales.  The  differences  stem  from  the  fact  that  d{x,y)  is  neither  a  Gaussian 
nor  is  d^x,  y)  equal  to  ^{x,  y). 

In  case  of  wavelets  that  are  equal  to  the  second  derivative  of  a  central  B-splines  [26,  8],  the 
two-dimensional  discrete  dyadic  wavelet  transform  decomposition  is  included  in  the 
mnltiscale  spline  derivative-based  one  from  our  previous  report.  To  realize  that,  let  us 
restate  the  filter  bank  implementation  of  the  multiscale  spline  derivative-based  transform 
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Figure  4:  Filter  bank  implementation  of  a  multiscale  spline  derivative-based  transform  with 
a  second  derivative  wavelet  for  m  G  [0,  M  —  1]:  (a)  Prefiltering,  (b)  postfiltering,  (c)  decom¬ 
position  and  (d)  reconstruction  modules. 


for  a  second  derivative  wavelet.  Figure  4  shows  the  building  blocks  of  a  filter  bank  with 
filters  specified  as 

GV)  =  eM^)(2jsin(f))^ 

F0(a;)  =  (cos  (f))^ 

T2{uj)  =  \H%u)\\  (32) 


and 


where  2p  — 1,  p  G  AT,  is  the  order  of  a  central  B-spline,  and  d  G  {1,2}. 

Note  that  the  decomposition  from  Figure  3  with  d—2  is  included  in  the  one  from  Figure  4. 


2.1.4  Enhancement  of  Mammographic  Features 

The  goal  of  the  method  in  this  section  is  to  adapt  specific  wavelet  based  enhancement 
schemes  for  distinct  mammographic  features,  and  then  combine  (fuse)  the  set  of  processed 
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Figure  5:  Overview  of  the  algorithm. 


images  into  an  enhanced  image.  The  mammographic  image  is  first  processed  for 
enhancement  of  microcalcifications,  masses,  and  stellate  lesions.  From  the  resulting 
enhanced  images,  the  final  enhanced  image  is  synthesized  by  means  of  image  fusion  [27]. 
Figure  5  presents  a  block  scheme  of  the  algorithm. 

Since  the  entire  processing  is  carried  out  within  the  wavelet  transform  framework,  both 
enhancement  and  fusion  are  implicit  (i.e.,  performed  in  the  wavelet  domain  only).  The 
algorithm  consists  of  two  major  steps:  (1)  wavelet  coefficients  are  modified  distinctly  for 
each  type  of  malignancy;  (2)  the  obtained  multiple  sets  of  wavelet  coefficients  are  fused 
into  a  single  set  from  which  the  reconstruction  is  computed.  The  devised  scheme  allows 
efficient  deployment  of  an  enhancement  strategy  appropriate  for  clinical  screening 
protocols:  enhancement  algorithm  is  first  developed  for  each  specific  type  of  feature 
independently,  and  the  results  are  then  combined  using  an  appropriate  fusion  strategy. 
(Note  that  it  is  possible  to  put  different  weights  on  features,  and  exclude  certain  features 
from  the  final  result.) 

Multiscale  spline  derivative- based  transform  with  a  second  derivative  of  a  central  B-spline 
wavelet  not  only  enables  the  inclusion  of  the  two-dimensional  discrete  dyadic  wavelet 
transform  decomposition  as  seen  in  Section  2.1.3,  but  also  allows  for  incorporation  of  a 
variety  of  methods  [7,  8,  28,  14,  9,  11,  12]. 

For  enhancement  of  microcalcifications,  second  derivatives  along  directions  of  x  and  y-axis 
are  added  to  form  an  approximation  to  a  Laplacian  of  Gaussian.  The  obtained  coefficients 
are  thresholded,  and  the  original  coefficients  at  the  corresponding  locations  multiplied  by  a 
gain  factor.  Similar  enhancement  through  detection  was  used  by  Strickland  and  Hahn 
[11,  12].  Within  the  multiscale  spline  derivative-based  transform  framework,  it  is  equally 
easy  to  obtain  voices  at  octaves  “2.5”  and  “3.5”  employed  in  [11,  12]:  central  B-spline 
properties  enable  computation  of  the  transform  at  any  integer  scale  [17]. 

Enhancement  of  circumscribed  masses  is  carried  out  by  applying  a  piecewise  linear 
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enhancement  function  [8]  to  wavelet  coefficients  at  scales  2^  through  2^  The  selection  of 
the  scale  range  was  based  upon  the  pixel  resolution  (116/rm)  of  digitized  mammograms  in 
the  University  of  Florida  database. 

Stellate  lesions  are  contrast  enhanced  according  to  the  observation  that  they  introduce  a 
distortion  into  a  radial  orientation  pattern  from  the  nipple  to  the  chest  wall  [28]. 
Approximately  steerable  [29]  second  derivative  wavelets  are  employed  in  conjunction  with 
their  Hilbert  transform  pairs  for  multiscale  orientation  analysis.  The  1-norm  of  differences 
between  local  orientation  and  average  orientation  within  a  sliding  window  is  used  as  an 
input  to  a  soft  thresholding  function  at  each  dyadic  scale  independently  [14]. 

Figure  6  demonstrates  the  results  of  processing  using  the  multiscale  analysis  contrast 
enhancement  algorithm  [9],  and  using  the  proposed  fusion  of  enhanced  features  method.  In 
this  example,  the  fusion  of  enhanced  features  emphasizes  the  appearance  of  a  mass  which  is 
surrounded  by  dense  parenchyma  of  the  breast.  Our  preliminary  results  suggest  that  this 
type  of  image  is  more  easily  interpreted  by  radiologists  compared  to  images  produced  via 
the  global  multiscale  analysis  contrast  enhancement  method  [9].  Note  also  that  the 
multiscale  analysis  based  global  contrast  enhancement  algorithm  [9]  can  readily  be 
incorporated  into  our  method  as  one  of  the  branches  before  the  fusion  module. 

The  choice  of  enhancement  parameters  controls  the  aggressiveness/subtleness  of  each 
resultant  enhancement  (i.e.,  the  prominence  of  the  targeted  features  with  respect  to  the 
surrounding  tissue).  The  structure  of  the  algorithm  is  also  well  suited  for  further 
refinements:  optimizations  can  be  performed  for  each  type  of  malignancy  alone,  and 
separately  for  the  fusion  module.  The  versatility  of  the  new  enhancement/fusion  paradigm 
can  provide  a  better  viewing  environment  to  facilitate  the  interpretation  of  mammograms. 

2.1.5  Implementation  Details 
Finite  Impulse  Response  Filters 

Since  all  two-dimensional  filters  used  in  the  filter  bank  implementations  of  the  transforms 
from  Section  2.1.3  are  x-y  separable,  we  will  begin  this  section  with  a  detailed  description 
of  FIR  filter  implementations  for  the  one-dimensional  discrete  dyadic  wavelet  transform 
implementation  from  Figure  1.  The  extension  to  two  dimensions  will  then  be 
straightforward . 

Let  us  refer  to  filters  applied  at  scale  2”®  as  filters  at  level  m-l-1,  and  let  filters  at  level  1 
(Equations  (10),  (14)  through  (16),  (27),  and  (32))  be  called  “original  filters,”  to 
distinguish  them  from  their  upsampled  versions.  As  an  input  to  the  filter  bank  from  Figure 
1,  we  consider  a  real  signal  s(n)  G  /^(Z),  n  G  [0,  iV  -  1].  Depending  on  the  length  of  each 
filter  impulse  response,  filtering  an  input  signal  may  be  computed  either  by  multiplying  the 
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(a)  (b)  (c) 


Figure  6:  (a)  Original  mammogram,  (b)  Contrast  enhancement  by  multiscale  analysis  (9). 
(c)  Enhancement  obtained  by  fusion  of  enhanced  features. 


33 


discrete  Fourier  transforms  of  the  two  sequences  or  by  circularly  convolving  s{n)  with  a 
filter’s  impulse  response.^  Of  course,  such  a  periodically  extended  signal  may  change 
abruptly  at  the  boundaries  causing  artifacts.  A  common  remedy  for  such  a  problem  is 
realized  by  constructing  a  mirror  extended  signal  [30] 


rs(-n-l)  ifne[-iV,-l] 

\  s{n)  if  n  G  [0,N  -  1], 


(33) 


where  we  chose  the  signal  Sme{'^)  fo  be  supported  in  [—N,N  1].  It  will  become  evident 

shortly,  that  mirror  extension  is  particularly  elegant  in  conjunction  with 
symmetric/antisymmetric  filters. 

Let  us  first  classify  symmetric/antisymmetric  real  even-length  signals  into  four  types  [31]: 


Type  I  /(n)  =  /(-n). 


Type  II  f{n)  =  f{-n  -  1), 


Type  III  f{n)  =  -/(-n). 

Type  IV  /(n)  =  -/(-n-l), 

where  n  E  [-N,N  -  1].  Note  that  for  Type  I  signals  the  values  at  /(O)  and  f{-N)  are 
unique,  and  that  for  Type  III  signals  the  values  at  /(O)  and  f{-N)  are  equal  to  zero. 

Using  properties  of  the  Fourier  transform,  it  is  easy  to  show  that  the  convolution  of 
symmetric/antisymmetric  real  signals  results  in  a  symmetric/antisymmetric  real  signal.  If 
a  symmetric/antisymmetric  real  signal  has  an  even  length,  then  there  always  exists  an 
integer  shift  such  that  the  shifted  signal  belongs  to  one  of  the  above  types. 

Now,  we  are  ready  to  examine  the  filter  bank  implementation  of  a  one-dimensional  discrete 
dyadic  wavelet  transform  from  Figure  1  with  filters  given  by  (10)  and  (14)  through  (16) 
driven  by  a  mirrored  signal  Sme{n}  at  the  input.  Let  the  number  of  levels  M  be  restricted 
by 

M  <  1  -f-  log2  ^  ■,  (34) 

^max  ^ 

where  L^ax  is  the  length  of  the  longest  original  FIR  filter  impulse  response. 

Each  FIR  filter  block  in  the  filter  bank  consists  of  a  filter  and  a  circular  shift  operator. 
Equation  (34)  guarantees  that  the  length  of  the  filter  impulse  response  does  not  exceed  the 
length  of  the  signal  at  any  block. 

Since  our  input  signal  Sme{n)  is  of  Type  II  and  noninteger  shifts  at  level  1  are  rounded  to 
the  nearest  integer,  it  follows  that  a  processed  signal  at  any  point  in  the  filter  bank  belongs 
®As  is  customary  in  image  processing,  we  use  circular,  rather  than  linear,  convolution. 
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to  one  of  the  types  defined  above.  This  means  that  filtering  a  signal  of  length  2N  can  be 
reduced  to  filtering  a  signal  of  approximately  one  half  of  its  length.  (For  Types  I  and  III, 

N  +  l  samples  are  needed.  However,  for  Type  III  one  needs  to  store  only  -  1  values 
because  zero  values  are  always  present  at  the  zeroth  and  (— A^)-th  sample  position). 
Implementation  is  particularly  simple  for  FIR  filters  designed  with  d  even  and  p  odd. 
Filters  are  of  Type  I  in  this  case,  so  the  signal  at  any  point  of  the  filter  bank  will  be  of 
Type  II.  An  FIR  filter  block  from  the  filter  bank  shown  in  Figure  1  can  therefore  be 
implemented  by 

L-1 

Fs,ni<^{n)  =  f{0)uji{n)  +  ^  /(^)[w//(n  —  2’^i)  +  un{n  +  2"^^)],  n  G  [0,  A'  —  1],  (35) 

i=l 

where 

(  u(-n-l)  ifnG[-Y,-l] 

Uii{n)  —  I  u{n)  ifnG[0,iV-l]  (36) 

[M(2A^-n-l)  ifnG[AA,f^], 

u{n)  is  an  input  signal  to  a  block,  f{n)  is  an  impulse  response  of  some  original  filter,  L  is 
the  length  of  the  filter,  and  N  is  the  length  of  an  input  signal  s(n)  to  the  filter  bank. 
Implementation  of  filters  bp{n)  used  for  prefiltering  and  postfiltering  represents  a  special 
case  of  (35)  with  m  =  0. 

A  filter  bank  with  the  above  implementation  of  blocks  and  signal  s(n)  at  the  input  yields 
equivalent  results  as  circular  convolution  for  s„e(^)  defined  by  (33).  In  addition  to 
requiring  one  half  the  amount  of  memory,  the  computational  savings  over  a  circular 
convolution  implementation  of  blocks  are,  depending  on  the  original  filter  length,  three  to 
four  times  fewer  multiplications  and  one  half  as  many  additions. 

A  similar  approach  can  be  used  for  other  filters.  However,  things  get  slightly  more 
complicated  in  this  case,  because  the  filters  are  not  of  the  same  type  and  the  signal 
components  within  the  filter  bank  are  of  distinct  types.  As  a  consequence,  an 
implementation  of  blocks  that  use  distinct  original  filters  may  not  be  the  same,  and  the 
implementation  of  blocks  at  level  1  may  differ  from  the  implementation  of  blocks  at  other 
levels  of  analysis. 

The  decomposition  blocks  at  level  1  can  be  implemented  by 

G-sfiu{n)  =  ^  g{i)[uji{n  —  z  —  1)  —  uii{n  +  z)],  n  G  [1,  A’  —  1], 
for  d  odd,  (35)  for  d  even, 

^-1 

2 

H-s,ou{n)  =  ^  h{i)[uii{n  —  z  —  1)  +  u//(n  +  z)],  n  G  [0,  N], 

1=0 
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for  p  even,  and  (35)  for  p  odd,  where  uji{l)  is  defined  by  (36),  g{n)  and  h{n)  are  impulse 
responses  of  the  filters  computed  from  (14)  and  (10),  respectively,  and  L  is  the  length  of 
the  corresponding  impulse  response. 

The  output  from  a  block  G-s(^)  at  level  1  is  of  Type  III  for  d  odd  and  of  Type  II  for  d 
even,  while  the  output  from  at  the  same  level  is  of  Type  I  for  p  even  and  of  Type  II 

for  p  odd. 

The  decomposition  blocks  at  subsequent  levels  m  G  [1,M-1]  can  be  implemented  by 
T-l 

G_,,„«(n)  =  Y,  s(t)Mn  -  2”(i  +  s))  -  «,(n  +  2"‘{i  +  s))],  n  £  |1,  iV  -  1], 

i=Q 

for  d  odd  and  p  even, 

^-1 

G_,,^u(n)  =  ^  g{i)[urj{n  -  2^(i  +  s))  -  u^n  +  2^(i  +  s))],  n  G  [0,  iV  -  1], 
i=0 

for  d  and  p  odd, 

L-1 

F^smujn)  —  /(O)uj(n)  +  -  2'^i)  +  ui{n  +  nG[0,  At],  (37) 

i—1 

with  /(n)  =  g{n)  for  d  and  p  even, 

H_,,mu{n)  =  ^  h{i)[ui{n  -  2^{i  +  s))  +  ur{n  +  2^{i  +  s))],  n  G  [0,  At],  (38) 

i=0 

for  p  even,  and  (35)  for  p  odd,  where 

(  u{-n)  if  n  G  [-y,  -1] 

ui{n)  =  <  uln)  ifnG[0,At]  (39) 

[  u{2N  —  n)  if  n  G  [N  +  1,  ^]. 

The  outputs  from  blocks  G'_s(2’"a;)  are  of  Type  III  for  d  odd  and  p  even,  of  Type  IV  for  d 

and  p  odd,  and  of  Type  I  for  d  and  p  even,  whereas  the  outputs  from  H^s{2^uj)  are  of 

Type  I  for  p  even  and  of  Type  II  for  p  odd. 

Next,  the  reconstruction  blocks  at  level  1  can  be  implemented  by 

L 

2 

Ksfiu{n)  =  k{i)[uiii{n  —  ^  +  1)  —  uiii{n  +  i)],  n  G  [0,  At  —  1], 

i=l 

for  d  odd,  (35)  for  d  even, 

L 

2 

Lsfiu{n)  =  '^l{i)[ui{n  -  z  +  1)  +  uj{n  +  z)],  n  G  [0,  At  -  1], 

irzl 
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for  p  even,  and  (35)  for  p  odd,  where 

'  -u{-n)  ifne[-f,-l] 

0  if  n  =  0 

uni{n)  =  \  uin)  ifnG[l,iV-l]  (40) 

0  if  n  = 

^  —u(2N  —  Ti)  if  n  G  [A^  +  1,  ^]) 

ui{n)  is  as  defined  by  (39)  and  k{n)  is  an  impulse  response  of  the  filter  from  (16).  Note 
that  both  outputs  from  blocks  Ks{u))  and  Ls{oj)  are  of  Type  II. 

The  reconstruction  blocks  at  subsequent  levels  can  be  implemented  by 

Ks,mu{n)  —  ^  k{i  +  l)[uiii{n  —  2'^{i  +  s))  —  uiji{n  +  2'”(i  +  s))],  n  G  [0,  N], 

i=0 

for  d  odd  and  p  even,  (37)  with  /(n)  =  k{n)  for  d  and  p  even, 

Ks,mu{n)  =  k{i  +  l)[u/y(n  —  2‘^{i  +  s))  —  uivin  +  2"*(f  +  s))],  n  G  [0,  Ai  —  1], 

2=0 

for  d  and  p  odd, 

for  p  even,  and  (35)  for  p  odd,  where  uiii{l)  is  given  by  (40), 

(  -u{-n-l)  ifnG[-Y,-l] 
ujv{n)  =  <  u{n)  if  n  G  [0, 77  -  1] 

[  -u{2N  -  n  -  1)  if  n  G  [Af, 

and  H^s,mu{n)  is  specified  by  (38).  We  observe  that  the  outputs  from  blocks  A's(2’^a;)  and 
Ls{2"^uj),  me[l,M  -  1],  are  of  Type  I  for  p  even,  and  of  Type  II  for  p  odd. 

When  we  compare  the  above  implementation  of  blocks  to  circular  convolution  driven  by  a 
mirrored  signal  s^e(n)  at  the  input,  we  observe  that  approximately  twofold  less  memory 
space,  three  to  four  times  fewer  multiplications  and  one  half  as  many  additions  are 
required.  (For  Type  I  signals  an  additional  sample  has  to  be  saved  because  two  values  are 
without  a  pair). 

Until  now,  we  have  talked  only  about  the  one-dimensional  case,  whose  filter  bank 
implementation  is  depicted  in  Figure  1.  Two-dimensional  transform  filter  bank 
implementations  (Figures  3  and  4)  are  not  only  comprised  of  x-y  separable  filters  solely, 
but  also  use  all  the  filters  from  Section  2.1.2.  Everything  presented  in  this  section  so  far  is 
therefore  directly  applicable  to  the  two-dimensional  case.  Filters  which  have  not  been 


37 


treated  yet  (i.e.,  ti{n),  t2(n),  and  bp{n))  can  all  be  realized  by  (35)  for  p  odd  or  m  =  0  and 
(37)  otherwise  (/(n)  denotes  an  impulse  response  of  any  of  the  above  mentioned  zero-phase 
filters  in  this  case). 

The  implementation  presented  in  this  section  performs  all  operations  in  the  spatial  domain, 
however,  one  could  also  implement  the  structures  shown  in  Figures  1,  3,  and  4  with  an 
input  signal  Sme{n)  (Equation  (33))  in  the  frequency  domain.  For  short  filter  impulse 
responses,  such  as  those  given  in  Tables  3,  4  and  5,  the  spatial  implementation  described  in 
this  section  is  certainly  more  efficient.  For  long  filter  impulse  responses,  however,  filtering 
is  faster  if  implemented  in  the  frequency  domain.  Additional  details  on  alternative 
implementation  strategies  can  be  found  in  [32]. 


Infinite  Impulse  Response  Filters 

Implementation  of  HR  filters  b~^{n)  which  were  introduced  in  Section  2.1.2  is  a  bit  more 
involved  than  the  one  encountered  in  the  previous  section.  Fortunately,  the  number  of 
different  cases  is  much  smaller  here:  possible  input  to  b~^{n)  in  filter  banks  from  Figures  1, 
3,  and  4  is  either  of  Type  II  or  of  Type  I.'^  We  will  use  ideas  and  a  few  results  from  [20]. 
Let  us  first  take  a  closer  look  at  the  system  function  B-^{z).  This  function  can  be  written 
as  a  cascade  of  terms 


E{z)  = 


1 


-a 


Z  —  -f  z~^ 

which  can  be  expressed  in  a  parallel  form  as 

-o;  /  1 


(1  —  q;2;”’-)(1  —  az) '' 


E{z)  =  - - ^  - - ZT  + 

1  —  \  1  —  az  ^ 


az 


(41) 


(42) 


where  a  and  ^  are  poles  of  the  causal  and  the  anticausal  filter,  respectively. 
The  impulse  response  of  this  term  can  be  written  as 


We  choose  to  implement  E{z)  in  a  cascade  form  and  therefore  extract  the  difference 
equations  from  (41): 


c^(n)  =  u{n)  +  ac'^{n  —  1)  n  =  1, 2, . . .  ,  N  —  1, 


(43) 


and 


c{n)  =  a  {c{n  +  1)  -  {n))  n  =  N -2,  N-3, . . .  ,0,  (44) 

"^Note  that  symmetry  types  in  this  section  slightly  differ  from  those  used  for  FIR  filters:  here,  mirror 
extended  signals  are  periodically  repeated,  so  that  they  stretch  from  -oo  to  oo. 
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where  u{n)  denotes  the  input  to  the  block,  c^{n)  is  the  output  from  the  causal  part,  and 
c{n)  stands  for  the  output  from  the  block. 

To  solve  (43)  and  (44)  we  need  boundary  conditions  c“^(0)  and  c(A^-l).  For  filters  b~^{n) 
in  filter  bank  implementations  from  Figures  1,  3,  and  4,  we  derive 


~  X/  ^  —  '^(0)  +  X/ 


a 


+  a 


2N-i 


io 


a 


2N 


-u{i)  ^  u(0)  +  ^  a'+^n(i),  (45) 


z~0 


and,  using  parallel  form  (42) 


a 


2N 


-uii)) 


—a 

1  - 


N-l 

(c+(iv-i)+  y] 

iz=N—l—io 


(46) 


where 


uii{n  mod  (2N))  if  n  >  0 

+  1)  iiiod  (2A^))  if  n  <  0, 


u{n)  if  n  e  [0,  -  1] 

u{2N-n-l)  if  n  G  [AT,  2A/' -  1], 


N  is  the  length  of  an  input  signal  to  the  filter  bank,  and  io<  N-l  is  selected  such  that 
falls  below  a  predefined  precision  threshold. 

For  orders  p  greater  than  three,  we  implement  Bp^{z)  as  a  cascade  of  terms  E{z)  with 
different  a’s.  Note  that  the  output  from  block  E{z)  is  always  of  the  same  type  as  the  input 
to  it. 


2.1.6  Summary 

We  extended  the  one-dimensional  discrete  dyadic  wavelet  transform  to  two  dimensions  and 
showed  that  in  case  of  wavelets  being  equal  to  the  second  derivative  of  a  central  B-spline 
function,  the  two-dimensional  discrete  dyadic  wavelet  transform  decomposition  is  included 
in  the  multiscale  spline  derivative-based  transform  decomposition.  This  attests  to  the 
fiexibility  of  the  latter  transform,  which  can  incorporate  a  variety  of  wavelet  based 
mammographic  image  enhancement  methods. 

The  multiscale  spline  derivative-based  transform  was  used  for  a  novel  image  enhancement 
scheme:  the  input  image  was  first  processed  for  enhancement  of  specific  mammographic 
features,  and  the  resultant  enhanced  images  were  fused  into  a  final  result.  The  described 
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scheme  represents  a  synthesis  between  our  early  global  methods  developed  for  an  overall 
improvement  of  image  quality  and  the  local  processing,  targeting  specific  mammographic 
features. 

The  employed  transform  is  highly  redundant,  and  efficient  implementation  is  of  particular 
importance.  We  presented  a  fast  filter  bank  implementation  that  takes  advantage  of  filter 
and  signal  symmetries.  Tests  of  our  method  included  convincing  quantitative  evidence  in 
terms  of  improved  visibility  and  detection  of  subtle  features. 
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2.2  Detection  of  Subtle  Masses  in  Mammography  via  Precision 
Multiscale  Analysis 
2.2.1  Introduction 

Breast  cancer  screening  by  mammography  is  the  most  effective  method  of  detecting  breast 
cancer  in  its  earliest,  most  treatable  stage.  On  average,  mammography  can  detect  tumors 
1.7  years  before  a  woman  can  feel  a  lump  via  self  examination.  Tumors  as  small  as  an 
eighth  of  an  inch  in  diameter  can  be  seen  with  mammograms  while  the  smallest  masses 
detectable  by  manual  exams  are  about  half  an  inch  across.  Early  detection  is  crucial  since 
survival  rate  has  a  roughly  inverse  relationship  with  the  stage  of  breast  cancer  at  detection, 
i.e.,  the  more  advanced  the  cancer  stage,  the  lower  the  survival  rate. 

Currently,  these  often  dense  and  complex  mammographic  images  are  examined  by 
radiologists  as  a  screening  modality.  However,  this  is  a  tedious  and  time-consuming  process 
since  only  one  in  a  thousand  images  will  show  features  that  are  of  possible  concern.  Today 
10  to  30%  percent  of  all  breast  cancers  are  missed  with  mammography,  as  currently 
practiced.  We  suggest  scale-space  search  algorithms  combined  with  the  emerging 
generation  of  digital  detector  technology  will  go  a  long  way  toward  reducing  false  negatives 
in  mammography  screening. 

Wavelet  algorithms  have  found  widespread  use  in  detection  applications  such  as  mass 
detection  in  mammograms  [33,  34,  35,  36,  37].  Wavelet  methods  are  preferred  over 
traditional  Fourier  analysis  because  the  local  nature  of  wavelets  allows  an  efficient 
representation  of  transients  by  a  small  number  of  components.  Traditional  Fourier 
transient  representations,  on  the  other  hand,  require  a  much  larger  number  coefficients  in 
order  to  prevent  sharp  discontinuities  in  the  spatial  domain  from  causing  ripple  artifacts  in 
the  frequency  domain.  Though  transients  are  typically  considered  to  be  a  phenomena  that 
occurs  in  time,  we  can  use  transient  formulation  to  represent  tumors  as  spatial  transients 
in  two-dimensional  mammograms.  Recent  studies  have  shown  that  indeed  mammographic 
features  are  well  represented  by  wavelet  coefficients  [9,  38]. 

In  wavelet  detection  theory,  Daubechies  first  pointed  out  that  the  wavelet  detector  should 
lie  in  the  same  spatial  neighborhood  and  in  the  same  frequency  band  as  the  feature  to  be 
detected  [39].  Therefore,  it  can  be  argued  that  the  dyadic  wavelet  transform  may  not  be 
appropriate  for  applications  where  subtle  detection  is  required.  That  is,  sampling  the  scale 
space  by  fixed  powers  of  two  may  not  provide  a  dense  enough  sampling  such  that  the 
wavelet  detector  lies  in  the  same  frequency  band  defined  by  the  feature.  In  fact,  in  this 
paper  we  show  that  the  dyadic  transforms  can  fail  in  mammography  applications.  In  the 
same  spirit,  Marco  proposed  to  use  an  M-band  wavelet  detector  to  obtain  a  denser  sampling 
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of  scale  space  [35].  In  this  paper,  we  describe  a  novel  translation-invariant,  continuous-scale 
wavelet  detector,  with  distinct  advantages  in  terms  of  automated  feature  detection. 
Previous  continuous  wavelet  algorithms  include  Rioul’s  [32]  and  Unser’s  methods  [40,  41]. 
In  Rioul’s  algorithm,  Q  new  mother  wavelets  are  defined  as  ss 
where  q  =  0, ...,  Q  -  1  and  ip^x)  =  h'^[n]4>{x  -  n)  in  each  octave.  An  octave-by-octave 

algorithm  is  applied  to  compute  “Q  voices  per  octave”  of  multiscale  analysis.  More  recently, 
Unser  developed  a  fast  algorithm  to  compute  the  wavelet  coefficients  for  Q  voices  within 
each  octave,  with  0{N)  complexity  per  voice.  However,  since  each  voice  corresponds  to 
similar  but  distinct  mother  wavelets,  it  is  not  straightforward  to  map  this  decomposition  to 
a  detection  algorithm.  In  terms  of  wavelet  detector  theory,  we  must  search  the  scale  space 
to  somehow  find  the  best  scale  for  a  signal.  Optimization  criterion  is  complicated  by  the 
fact  that  any  change  in  scale  might  require  a  change  in  the  mother  wavelet.  In  this  paper, 
we  describe  an  algorithm  to  calculate  the  wavelet  coefficients  for  the  same  mother  wavelet. 
In  our  method,  a  continuous  scale  wavelet  detector  is  incorporated  into  a  Constant 
False-Alarm  Rate  (CFAR)  detector  to  improve  performance.  A  CFAR  detector  utilizes  the 
statistics  of  a  target  and  its  background  to  separate  the  target  from  its  background  and  is 
widely  used  in  Automatic  Target  Recognition  applications  [42].  We  show  that  such  a 
detector  is  able  to  effectively  separate  a  mass  from  its  complex  background  using  the 
respective  statistics  of  each  region. 

This  chapter  is  organized  as  follows.  In  section  2,  we  derive  a  arbitrary  scale  discrete 
wavelet  decomposition  algorithm.  Next,  in  section  3,  we  describe  a  arbitrary  scale  wavelet 
CFAR  detector.  In  section  4,  we  show  that  the  best  scale  to  detect  a  tumor  is  related  to 
the  size  of  each  distinct  tumor.  Then  we  show  the  results  of  applying  our  algorithm  to 
difficult  mammography  cases  of  breast  cancer.  The  ability  to  find  and  process  features  at 
an  arbitrary  scale  is  shown  to  have  significant  advantages  over  standard  dyadic  algorithms. 
Finally,  in  section  5,  we  summarize  our  results. 

2.2.2  A  Continuous  Scale  Discrete  Wavelet  Transform 

A  wavelet  transform  is  a  linear  operation  that  projects  a  signal  onto  a  set  of  basis 
functions,  which  are  dilated  versions  of  a  mother  wavelet.  The  mother  wavelet  is  a  function 
Ip  e  that  satisfies  <  +oo-  This  requires  that  the  wavelet  has  sufficient 

decay,  and  that  ip{x)dx  =  0.  The  continuous  wavelet  transform  of  a  function  at  scale  a 
and  shifting  parameter  b  is  defined  as 

W/(a,  b)^  J  f{x)i;  dx. 
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Figure  7:  Data  acquisition  and  initial  condition:  (a)  Block  diagram  of  digitization  a  con¬ 
tinuous  signal  f(x)  is  sampled  after  a  low-pass  operation,  (b)  initial  condition  for  subsequent 
wavelet  decomposition. 

Of  course,  it  has  been  shown  that  a  wavelet  transform  satisfies  energy  conservation  and  can 
be  perfectly  reconstructed  from  its  wavelet  transform  coefficients  [16,  43].  The  Arbitrary 
Scale  Discrete  Wavelet  Transform  (ASDWT)  is  defined  as  the  CWT  sampled  along  the 
shifting  parameter  b, 

W/(a,n)  =  J"  /(x)^ 

In  practice,  data  is  spatially  sampled  and  quantized.  Furthermore,  we  assume  that  a  signal 
f{x)  is  sampled  at  its  Nyquist  rate  fs  so  that  it  may  be  recovered  from  its  discrete  samples. 
Since  f{x)  is  band-limited,  W/(a,  6)  is  also  band-limited,  and  can  be  recovered  from  its 
samples  W/(a,  n). 

Scaling  Space 

The  finite  resolution  of  the  sampled  signal  g[n]  puts  a  limit  on  the  finest  scale  that  can  be 
computed.  Let  us  now  introduce  a  compactly  supported  scaling  function  (l){x)  that  satisfies 
the  following  conditions  [40]  : 

1.  Riesz  basis  condition:  A  <  +  2A:7r)p  <  B,  where  A  and  B  are  strictly 

positive  constants; 

2.  Order  property:  $(0)  =  1,  <F^”*^(2A:7r)  =  0,k  £  Z,k  ^  0,  for  m  =  0, ...,  —  1; 

3.  Two-scale  relation:  (f)  (|)  =  Ylkez  h[k](l){x  —  k)] 

4.  Frequency  property:  $(a;)  7^  0,  for  Iwj  <  tt,  and  $(w)  should  be  very  small  for 

IcjI  >  TT. 

Constraint  1  implies  that  <p{x  —  k)  is  a  Riesz  basis  of  the  subspace  defined  by 
'^<p  =  {f{x)  =  '^c[k]4>{x  -  k)  CGI2}  foTf{x)e'V^ 

kez 
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and  that  is  a  well  defined  (closed)  subspace  of  £2  [44].  Property  2  implies  that  (p{x) 
reproduces  all  polynomials  of  degree  N-1  [45].  Constraint  3  is  the  well-known  two-scale 
relation.  Condition  4  shows  that  (j){x)  can  be  an  anti-aliasing  filter. 


Initial  Condition 

In  the  traditional  DWT  analysis  literature,  it  is  generally  assumed  that 


/OO 

f{x)4>{n  -  x)dx. 

-00 


(48) 


However,  in  practice,  g[n]  is  the  output  of  a  physical  acquisition  device  shown 
schematically  in  Fig  7  (a).  The  impulse  response  of  an  acquisition  device  is  ^(^x),  where 
^(a;)  is  typically  a  lowpass  filter.  The  approximation  of  ^{x)  in  space  is 

^{x)  ^  i{x)  =  '^r[k]^(x  -  k),  (49) 

k 

where  r[k]  is  the  projection  of  ^{x)  onto  the  scaling  space.  Let  g'[n\  =  f{x)(f){n  —  x)dx. 

The  relation  between  g'[n]  and  g[n\  is 


9[n] 


/OO 

f{x)^{n  —  x)dx 

-OO 

/OO 

r[k]4>{n  —  X  —  k))dx 
k 

/OO 

f{x)(f>{n  —  X  —  k)dx) 

K 

=  Y^r[k]g'[k], 


(50) 


Therefore,  g'[n]  satisfies  the  traditional  initial  condition  and  can  be  calculated  from  g[n\  by 

9'[n]  =  Yl9[k]'r'~^{n  -  k).  (51) 

k 

It  is  commonly  assumed  that  ^{x)  is  an  ideal  lowpass  filter.  An  ideal  lowpass  filter  may  be 
approximated  as  a  cardinal  spline  filter  [46].  For  example,  the  scaling  function  is  a  spline  of 
order  n  —  1, 


(52) 


Therefore,  r~'^[k]  =  b^[k].  Notice  that  Equation  (51)  is  similar  to  projecting  a  signal  onto  a 
scale  space  as  described  in  [40]  [41],  although  the  theoretical  justification  for  our  case  is 
quite  different. 
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ASDWT 

A  discrete  approximation  of  f{x)  at  scale  a  is  defined  as 


Sa/N  =  J  fi^)^ 

It  can  be  shown  from  the  two-scale  relation  and  Equation  (53)  that 

;F{Sj„/[nl}  =  H(Vw)T{Si,l[n]}-  (54) 

From  the  initial  condition,  Si/[/i]  =  where  a  =  1  is  the  finest  scale  of  analysis.  As  a 
increases,  Saf[n]  becomes  a  coarser  representation  of  f{x).  When  a  goes  to  infinity,  Sa/[n] 
approaches  a  constant  equal  to  the  average  value  of  f{x).  In  practice,  the  length  of  g[n] 
defines  the  coarsest  scale. 

When  $(cc;)  /  0  for  a;  G  [-7r,7r],  it  is  straightforward  to  show  that  there  is  one  and  only  one 
f{x)  such  that 

g'^ix)  =  f[x)  *  (i){x).  (55) 


Theorem  1  Given  g[n],  for  a  >  sq,  if^ai}^)  *5  band-limited  between  [-7r,7r],  the  ASDWT 
can  be  calculated  by 


where 


and 


Waf{a,  n)  =  XI ’ 

kez 


P%a,Lv) 


^(o;) 


P'^{a,u;) 


(56) 


(57) 


Proof:  Note  that  'ipa{x)  and  p'^{a,x)  are  band-limited  at  [— tt,  tt].  And 
p{a,n)  =  p^{a,x)\x=n-  From  Equation  (57), 

^a{co)  =  P‘'{a,ujMio)  (58) 


'^a{x) 


u)du. 


(69) 


45 


If  we  substitute  Equation  (59)  into  the  previous  ASDWT  definition  (Equation  (47)),  we 
obtain 


Therefore,  we  can  claim  that  Equation  (61)  is  indeed  true.  Substituting  Equation  (61)  into 
Equation  (60),  we  obtain 

Waf{a,  n)  =  ^p(a,  k)g{n  -  k). 
kez 


Based  on  Theorem  1,  p{a,k)  can  be  calculated  by  sampling  IFFT  as  long  as 

^o(i^)  is  band-limited  between  [-7r,7r].  Once  p{a,  k)  is  available,  the  sample  of  CWT  of  a 
signal  can  be  calculated  [32]. 

Note  that  the  mother  wavelet  must  be  band-limited  so  that  we  can  obtain  wavelet 
coefficients  by  discrete  wavelet  transform  using  the  same  mother  wavelet.  Note  that  a 
wavelet  ip{x)  is  not  required  to  be  compact  in  the  frequency  domain.  However,  if  4f(a;)  is 
local,  then  we  can  select  an  sq  such  that  ^so(w)  =  has  an  acceptable 
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approximsition  error  onto  the  scaling  space.  The  approximation  error  is  defined  as  the  ratio 
of  the  energy  of  $(aic)  ontside  [-x,  x]  to  the  total  energy.  The  ratio  decreases  as 

increases. 

Next  we  show  that  our  4.(x)  is  also  a  mother  wavelet.  Since  satisfies  the 
“admissibility  condition,” 


/OO 

-OO 


OO. 


Substituting  co 


with  sou;  and  multiplying  both  side  of  the  equation  with  Rect{u^),  we  have 


“^|§(sow)pda;  < 


OO 

OO. 


(62) 


Therefore,  satisfies  the  “admissibility  condition.”  Wavelet  coefficients  at  arbitrary 

scale  a  =  ssq  can  be  calculated  directly  from  Equation  (56). 

When  little  is  known  about  the  best  scale  to  detect  a  certain  feature,  it  is  desirable  to 
search  the  scale  space  at  a  coarse  discretization  (perhaps  dyadically)  before  carrying  out  a 
finer  search.  Once  an  so  is  selected,  it  may  be  considered  as  a  “voice”  and  F(2  so,u;)  can 

be  calculated  as  [32] 


L-l 


P(2^so,a;)  =  P{so,2^uj)(Y[H{2>^uj)). 

k^O 


(63) 


Thus,  Eq(54),  (56)  and  (63)  yield 

^{Wsiso/}  =  P{so,2^oj)P{S2if}. 


(64) 


In  this  study,  the  wavelet  is  a  second  derivative  of  a  spline  of  order  5  (n  4,  d  2).  The 
scaling  function  and  mother  wavelet  are  given  by 


$(w) 


sm(|) 

2 


n 


4r(a;)  =  {juj^ 


n-{-d 


and 


H(u,)  =  e«"cos”  (I)  ,  P(»o.w)  = 


^d2n+2d 

(I) 


xo  =  l  when  n  is  odd  and  level  L  =  0.  Otherwise,  xq  =  0.  xi  =  |  when  d  is  odd  and  level 
L  =  0,  and  Xi  =  0. 

In  next  section,  we  will  apply  this  ASDWT  in  tumor  detection  in  difficult  mammography 
cases  of  breast  cancer. 
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2.2.3  A  CFAR  Wavelet  Detector 

In  this  section,  we  first  describe  a  novel  continuous  scale  wavelet  detector.  Then  we 
introduce  the  CFAR  detection  theory  and  show  that  our  wavelet  detector  is  also  a  CFAR 
detector  under  certain  assumptions.  Later,  we  show  the  results  of  applying  our  algorithm 
to  difficult  mammography  cases  of  breast  cancer.  The  ability  to  find  and  process  features 
at  an  arbitrary  scale  is  shown  to  have  significant  advantages  over  standard  dyadic 
algorithms.  Finally,  we  summarize  our  results. 

Wavelet  Transient  Detector 

We  pose  mass  detection  in  mammograms  as  a  transient  signal  detection  problem.  Let’s 
consider  hypotheses  of  the  form 

Ho  :  y{x)  =  n{x), 

Hi  :  y{x)  =  f{x)+n{x). 

Where  x  is  the  location,  y{x)  is  the  observed  signal,  f{x)  is  the  underlying  signal  to  be 
detected,  and  n{x)  is  additive  noise.  Based  on  the  observation  y,  it  must  be  decided 
whether  the  tumor  is  present  (ifi)  or  not  (Hq). 

We  assume  that  a  signal  f{x)  is  essentially  spatial-localized  and  band-limited  ^  (to 
spatial-frequency  area  [xj.^Xf^]  x  {[-cof^,  -w/J  U  [cof^,u!f^]}),  as  shown  in  Figures.  Then  we 
project  f{x)  onto  a  wavelet  basis  at  scale  s^: 

^sif{b)  =  j  f{x)tl^siib-  x)dx.  (65) 

The  advantage  of  applying  the  wavelet  representation  to  this  detection  problem  is  a  result 

of  the  well  known  fact  that  there  are  only  a  few  significantly-different-from-zero  wavelet 

coefficients  for  an  essentially  spatial-localized  and  band-limited  signal  [39].  The 

spatial-frequency  area  of  a  wavelet  V’si  is  shown  in  Figure  8.  As  scale  s  changes,  the 

spatial-frequency  area  of  '0s  moves  vertically.  As  the  shifting  parameter  h  in  Equation  (65) 

changes,  the  spatial-frequency  area  of  0s  shifts  horizontally.  When  the  spatial-frequency 

area  of  0s(^  -  x)  overlaps  with  that  of  f{x),  wavelet  coefficients  Wsf{b)  may  be 

significantly  different  from  zero.  Note  that  when  we  project  the  observed  signal  y{x)  on 

wavelet  base  0s. ,  we  filter  out  the  noise  that  lies  outside  the  passband  of  05^  and  increase 

the  SNR.  Therefore,  we  can  consider  the  process  of  decomposing  an  observed  signal  y(x) 

onto  the  best  scale  s*  as  projecting  that  signal  onto  a  basis  to  maximize  the  SNR. 

^We  want  to  detect  a  mass  whose  size  is  in  a  certain  range.  In  later  discussion,  we  show  that  the  size  of  a 
mass  is  corresponding  to  a  certain  scale  [frequency].  Therefore,  it  is  reasonable  to  assume  a  profile  of  a  mass 
is  spatial-localized  and  band-limited. 
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Mallat  has  experimentally  found  that  the  wavelet  coefficients  can  be  modeled  with  a 
two-parameter  generalized  Laplacian  distribution  [47]: 

p(m)  =  (66) 


where  the  constant  C  is  adjusted  in  order  to  have  p[u)du  —  1.  The  coefficients  a  and 
(3  are  directly  related  to  the  first  and  second  moment  of  the  wavelet  coefficients,  defined  as 

/+00  r-\-oo 

\u\p{u)du  and  m2  =  /  u^p{u)du. 

-00  — 00 


We  can  compute  a  and  (3  using 


13  = 


J 


_  m2T{l/l3) 

"  r(3//?)  ’ 

where 

fte)  -  r(2/x)^ 

r(3/x)r(i/x)’ 

roo 

r(x)  =  /  u^-^e-^du. 

Jo 

Therefore,  the  pdf  of  wavelet  coefficients  under  hypothesis  Ho  is  Po(m)  =  .  The 

pdf  of  wavelet  coefficients  under  hypothesis  Hi  is  Pi(m)  =  \  and 

m  =  MAX{WsJ).  The  “Likelihood  Ratio  Test”(LRT)  becomes 

A{u)  =  log{po{u))  -  log{pi{u))  (®’^) 
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As  in  Figure  9,  when  A(u)  <  X,  u  E  Rq  and  u  <T.  When  A{u)  >  X,  u  ^  R\  and  u  >  T. 
Note  that  no  matter  where  we  set  the  threshold  T,  we  will  occasionally  make  a  wrong 
decision.  When  Hq  is  true  and  we  choose  Hi,  we  call  it  “error  of  the  first  kind.”  Its 
probability  Qo  is 


it  is  represented  by  the  area  under  the  curve  of  po  to  the  right  of  T  (the  shadowed  area  in 
Figure  9).  We  also  call  it  a  “False  alarm.”  On  the  other  hand,  we  might  choose  Hq  when 
Hi  is  true.  That  is  an  “error  of  the  second  kind.”  Its  probability  Qi  is 

Qi  =  /  pi{u)du. 

Of  course,  the  value  of  T  depends  on  how  much  these  mistakes  cost.  It  also  depends  on  the 
prior  probability  of  Hi.  For  example,  if  we  know  that  the  family  of  a  woman  has  a  history 
of  breast  cancer,  she  is  at  high  risk  of  getting  breast  cancer  and  we  should  lower  the 
detection  threshold  T. 

After  the  hypothesis  testing,  we  can  utilize  the  geometric  information  to  further  eliminate 
false  alarms.  Our  ASDWT  detector  consists  of  the  following  procedures: 

1.  Find  the  best  scale  by  systematically  searching  the  scale  space,  from  coarse  to  fine. 
(The  definition  of  “best  scale”  is  given  later.) 

2.  Choose  the  maximum  ASDWT  coefficients  from  the  best  scale. 

3.  Eliminate  the  wavelet  maxima  which  are  less  than  some  threshold.  The  value  of  the 
threshold  is  determined  by  the  desired  false  alarm  rate. 

4.  If  local  curvature  analysis  of  the  remaining  targets  indicates  a  symmetric  bump,  the 
center  of  the  bump  is  considered  to  be  the  site  of  a  potential  mass  (see  discussion 
below). 

Each  pixel  that  is  larger  than  a  threshold  (determined  by  the  false  alarm  rate)  is  a 
candidate  target.  Local  curvature  analysis  seeks  to  verify  that  each  candidate  target  is  a 
symmetric  bump  with  a  single  local  maximum.  Asymmetric  bumps  are  artifacts  due  to 
edge  effects  or  extended  linear  features  in  the  dense  tissue.  The  characteristic  star  patterns 
in  the  figures  are  created  using  a  simple  algorithm  that  compares  each  pixel  with  it 
neighboring  pixels  and  decides  if  the  pixel  is  a  local  maximum  in  any  direction. 
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(a) 


Figure  10:  A  conventional  CFAR  detector:  (a)  Standard  interior  and  exterior  regions,  (b) 
3-D  plot  of  the  equivalent  filter. 


A  CFAR  Detector 

In  detection  theory,  when  the  prior  probabilities  of  the  two  hypotheses  Hi  and  H2  and  the 
costs  associated  with  each  type  of  error  are  known,  a  threshold  T  should  be  chosen  to 
minimize  the  average  risk  by  using  the  “Bayes  criterion.”  If  the  prior  probabilities  are 
unknown,  the  “minimax  criterion”  should  be  used  to  select  the  threshold  T.  When  not 
only  the  prior  probabilities  but  also  the  costs  are  difficult  to  estimate,  one  should  select  an 
acceptable  false  alarm  rate  and  try  to  maximize  the  detection  rate.  That  is  the 
“Neyman-Pearson  criterion.”  Goldstein  proposed  a  automatic  detection  procedure  which 
maintained  a  Constant  False- Alarm  Rate  (CFAR)  [48].  He  showed  that  the  false-alarm  rate 
depended  on  a  detection  threshold  and  the  number  of  data  used  to  estimate  statistic 
parameter.  Based  on  a  clutter  distribution  assumption  and  test  statistic  theory,  he 
provided  a  procedure  to  adjust  the  detection  threshold  adaptively  in  order  to  keep  the  false 
alarm  rate  constant. 

The  CFAR  detection  algorithm  is  widely  used  in  military  Automatic  Target  Recognition 
(ATR)  applications.  A  one-parameter  CFAR  detector  is  a  cell-averaging  detector  controls 
the  detection  threshold  for  a  specific  cell  (pixel)  based  on  the  estimation  of  a  sufficient 
statistic  of  the  clutter  [42].  It  works  well  under  the  assumption  that  the  background  noise  is 
statistically  homogeneous.  During  hypothesis  testing,  a  test  pixel  Xt  is  set  to  be  the  center 
of  test  region  Rt-  The  mean  of  test  region  Rf  is  calculated.  Then  the  clutter  (background) 
intensity  mean  is  estimated  in  a  local  region  Rc  distantly  surrounding  pixel  Xt,  which  is 
shown  in  Figure  10  (a).  Their  difference  is  compared  to  a  threshold  Tcfar  to  determine 
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whether  pixel  Xt  is  a  target.  This  one-parameter  CFAR  detection  can  be  formulated  as 

*  {i,j)eRc 

And 

^target  rp 

y  ^non— target  ^CFAR^ 

where  x{i,j)  is  the  intensity  of  pixel  at  location  Nt  and  Nc  are  the  number  of  pixels 
in  region  Rt  and  Re-  Note  that  only  the  mean  information  is  utilized  m  this  one-parameter 
CFAR  detector.  In  order  to  drop  the  homogeneous  requirement,  we  can  adapt  a  threshold 
to  a  local  environment.  A  two-parameter  CFAR  detector  utilizes  both  first  and  second 
order  statistics  in  hypothesis  testing,  and  is  given  by 

t‘':  =  W  T, 


O-c  = 


Y 

{i,j)ERc 


•> 


y  =  Xt-  Me, 


y  ^nm-\arget  CTcTcFAR- 

From  a  signal  processing  perspective,  outputs  y  of  a  CFAR  detector  can  be  obtained  by 
filtering  the  image  with  a  kernel  shown  in  Figure  10  (b).  The  abrupt  changes  in  this  kernel 
will  cause  undesirable  ripples  outside  its  main  lobe  in  the  frequency  domain.  We  propose  to 
use  a  wavelet  which  is  a  second  derivative  of  a  spline  of  order  5  (n=4,  d=2).  As  shown  in 
Figure  11,  this  ASDWT  detector  is  a  weighted  CFAR  detector.  When  calculating  the 
statistics  of  the  test  area,  the  pixels  in  the  center  of  the  testing  area  are  weighed  more 
heavily.  While  estimating  the  statistic  of  the  surrounding  area,  pixels  both  near  the  test 
area  and  far  away  from  the  test  area  are  weighed  less  heavily.  While  the  shape  of  our 
operator  resembles  the  classical  center  surround  operators,  it  should  be  remembered  that 
image  statistics  are  gathered  in  the  positive  and  negative  valued  regions  independently. 
Therefore,  the  proposed  detector  combines  a  wavelet  formulation  with  the  classical  theory 

of  CFAR  detectors. 

Next,  we  will  apply  our  wavelet  CFAR  detector  for  mammographic  tumor  detection. 
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Results  of  Mammographic  CFAR  Detector 

In  this  section,  we  use  a  digital  radiograph  of  a  phantom  to  demonstrate  that  continuous 
scale  is  necessary.  We  show  that  our  pd/ assumption  is  reasonable.  And  we  apply  our 
detector  to  actual  mammograms  obtained  from  Shands  Hospital  in  Gainesville,  FL. 

X-ray  images  of  an  RMI  model  156  phantom  (Radiation  Measurements  Inc.,  Middleton, 

WI)  were  investigated.  The  RMI  phantom  contains  five  masses  and  six  fibrils  which  mimic 
two  objects  of  interest  in  mammography.  The  radiographic  image  quality  of  this  phantom, 
as  well  as  the  total  number  of  objects  which  are  deemed  to  be  visible,  form  a  key 
component  of  the  mammography  accreditation  program  administered  by  the  American 
College  of  Radiology  (ACR).  The  RMI  phantom  is  designed  so  that  at  least  one  object  (i.e. 
smallest  mass)  in  each  category  is  not  visible,  and  one  object  is  at  the  borderline  of  the 
visibility  threshold  [49]. 

Radiographic  images  of  the  RMI  phantom  were  obtained  using  a  Digital  Spot 
Mammography  (DSM)  system  attached  to  a  Lorad  Breast  Biopsy  System  (Lorad 
Corporation,  Danbury,  CT).  Due  to  the  limit  of  the  detector  area,  only  regions  containing 
masses  at  the  right-lower  corner  were  captured  (with  technique  factors  of  22  kVp  and  16 
mAs)  as  shown  in  Figure  12(a).  The  smallest  mass  in  the  schematic  representation  of  the 
insert  of  the  phantom  (Figure  12(e))  can  not  be  seen  using  conventional  window  and 
leveling  at  12-bit  contrast  resolution. 

First,  a  wavelet  detector  at  scale  64  was  applied.  A  threshold  was  chosen  such  that  there 
are  no  false  alarms.  As  shown  in  Figure  12  (c),  both  the  largest  mass  and  the  medium-size 
mass  were  detected  while  the  smallest  mass  was  missed.  Then  we  reduced  the  scale  of  the 
detector  to  scale  46.4  and  found  that  all  the  masses  were  detected  with  zero  false  alarms  as 
shown  in  Figure  12  (e).  When  the  scale  is  set  to  a  dyadic  power  32,  we  cannot  detect  the 
smallest  mass  without  a  false  alarm  as  shown  in  Figure  12  (f).  What  is  more,  the  detector 
cannot  even  detect  the  center  of  the  largest  mass  correctly.  When  we  look  at  the  wavelet 
coefficients  shown  in  Figure  12  (c),  (d)  and  (e),  the  smallest  mass  can  be  seen  at  scale  46.4 
and  is  clearly  less  evident  at  scale  32  and  64.  Since  the  dyadic  wavelet  detector  failed  to 
detect  all  masses  without  false  alarms,  a  more  arbitrary  and  flexible  scale  of  analysis  is 

necessary. 

The  non-dyadic  scale  46.4  in  the  above  paragraph  was  found  by  using  an  estimated  range 
of  mass  sizes  to  bracket  a  search  for  an  effective  scale.  An  effective  scale  provides  the 
characteristic  star  pattern  in  the  output.  We  will  further  discuss  the  optimal  scale  search  m 

next  chapter. 

Having  proved  the  case  using  phantoms,  we  now  apply  our  detector  to  some  real 
mammograms  with  dense  tissue  obtained  from  Shands  Hospital  in  Gainesville,  FL.  In  this 
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Figure  12:  Continuous  scale  is  important,  (a)  Digital  radiograph  of  the  RMI156  phantom, 
(b)  Schematic  representation  of  mammographic  features  within  the  phantom,  (c)  Detector 
output  obtained  at  scale  64,  super-imposed  on  the  phantom  image,  (d)  Wavelet  coefficients 
at  scale  64.  (e)  Detector  output  obtained  at  scale  46.4.  (f)  Wavelet  coefficients  at  scale  46.4. 
(g)  Detector  output  obtained  at  scale  32.  (h)  Wavelet  coefficients  at  scale  32. 


Case 

Shape 

Margins 

Density 

Size 

Biopsy 

Difficulty 

lcc046 

Irregular 

Indistinct 

High 

10mm 

Malignancy 

Mod.  Difficult 

lml015 

Round 

Obscured 

High 

15mm 

Malignancy 

Difficult 

lcc002 

Irregular 

Spiculated 

High 

10mm 

Malignancy 

Difficult 

lcc004 

Irregular 

Spiculated 

High 

15mm 

Malignancy 

Difficult 

Table  7:  Description  of  mammograms. 


preliminary  study,  we  tested  our  algorithm  using  cases  that  a  mammography  expert 
regarded  as  likely  to  be  missed  by  radiologists.  These  mammograms  contain  lesions  of 
varying  density  as  described  in  Table  7.  In  figure  13,  we  show  histograms  of  wavelet 
representations  of  mammographic  background  in  case  lcc0046  and  lml015  and  histograms 
of  the  corresponding  tumor  wavelet  representations  .  All  those  histograms  are 
exponentially  decay  as  what  we  assume.  Because  we  have  a  limited  number  of  pixels,  the 
histograms  are  noisy.  In  Figure  14,  15,  16  and  17,  we  show  some  initial  results.  All  masses 
(biopsy  proven)  were  correctly  detected.  The  star  patterns  in  each  figure  are  added  by  the 
local  curvature  analysis  discussed  earlier.  All  symmetric  star  pattern  indicate  possible 
masses  that  turned  out  to  be  actual  masses  as  verified  by  biopsy  procedures. 

In  the  next  section,  we  will  discuss  how  to  select  the  best  scale  for  tumor  detection. 

2.2.4  Best  Scale  Selection 

In  this  section,  we  describe  a  method  to  find  the  optimal  scale  to  segment  a  potential  mass 
from  its  background,  using  a  continuous  multiscale  analysis.  First,  a  suspicious  object  is 
identified,  which  may  or  may  not  contain  a  lesion.  Next,  we  search  the  scale  space  to  find 
the  scale  at  which  the  contrast  of  the  object  is  maximum.  Finally  we  will  apply  this 
scheme  to  some  real  mammograms. 

In  the  next  subsection,  we  will  discuss  how  to  identify  1-D  image  objects. 

Relationships  between  Scale  and  Objects 

Singularities  and  irregular  structure  often  carry  the  most  important  information  about  an 
object.  For  example,  the  discontinuities  in  the  intensity  provide  the  locations  of  the  object 
contours,  which  are  particularly  meaningful  for  recognition  purposes.  The  local  regularity 
of  images  can  be  characterized  by  wavelet  representations.  In  this  section,  we  will  first 
introduce  the  Lipschitz  exponents,  which  is  a  measure  of  the  local  regularity  in 
mathematics.  Then  we  will  review  the  Lipschitz  exponents  and  their  characterization  with 
the  wavelet  transform.  Finally,  we  will  explain  how  to  locate  an  object  to  find  the  best 
scale  to  detect  that  object. 
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(b) 

Figure  13:  Normalized  histograms  of  tumor  wavelet  representation  (dotted  line)  and  mam- 
mographic  background  wavelet  representation  (solid  line),  (a)  Mammogram  case  lcc046.  (b) 
Mammogram  case  lml015. 


(c) 


Figure  14:  Detection  result:  (a)  Mammogram  case  lcc046,  (b)  the  maxima  of  outputs  of  the 
wavelet  CFAR  detector,  (c)  mass  enhanced  by  window  and  leveling. 
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Figure  16:  Detection  result:  (a)  Mammogram  case  lml015,  (b)  the  maxima  of  outputs  of  the 
wavelet  CFAR  detector,  (c)  mass  enhanced  by  window  and  leveling. 
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(c) 


Figure  17:  Detection  result:  (a)  Mammogram  case  lcc002,  (b)  the  maxima  of  outputs  of  the 
wavelet  GEAR  detector,  (c)  mass  enhanced  by  window  and  leveling. 
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In  matliGniatics.  the  local  regularity  is  often  measured  with  Lipschitz  exponents,  which  is 
defined  as  [50,  15]; 

•  Let  n  be  a  positive  integer  and  n  <  a.  <  n-\-l.  A  function  f{x)  is  said  to  be  Lipschitz 
a,  at  xq,  if  and  only  if  there  exists  two  constants  K  and  ho,  and  a  polynomial  of 
order  n,  Pn{x),  such  that  for  h  <  ho 

|/(a:o  +  h)  —  Pn{h)\  <  K\h\  ■ 

•  Let  0  <  a  <  1.  A  function  f{x)  is  uniformly  Lipschitz  a  over  an  interval  ]a,b[  if  and 
only  if  there  exists  a  constant  K  such  that  for  any  {xo,Xi)  ^]a,b[^ 

\f{xo)  -  /(a:i)|  <  K\xo-xi\^. 

•  We  call  Lipschitz  regularity  of  f{x)  and  xo,  the  superior  bound  of  all  values  a  such 
that  f{x)  is  Lipschitz  a  at  xq. 

•  We  say  the  a  function  is  singular  at  Xo,  if  it  is  not  Lipschitz  1  at  xq. 

If  /(x)  is  discontinuous  but  bounded  in  the  neighborhood  of  Xq,  its  uniform  Lipschitz 
a  =  1.  If  /(x)  is  differentiable  at  xq,  then  it  is  Lipschitz  a  =  1.  The  larger  uniform 
Lipschitz  a  is,  the  more  “regular”  /(x)  at  xq  is.  We  can  measure  the  Lipschitz  exponents 
from  the  evolution  across  scales  of  the  absolute  value  of  the  wavelet  transform.  We  suppose 
a  wavelet  ■0(x)  is  continuously  differentiable  and  decays  asymptotically  as  Let 

0  <  Of  <  1.  A  function  is  uniformly  Lipschitz  a  over  ]a,  6[,  the  wavelet  transform  satisfies 
[50]: 

toS2|W./(3:)|  <  log2lK)  +  alog2(s).  (68) 

Note  that  a  discrete  signal  is  of  finite  resolution,  which  we  normalize  to  1.  We  cannot 
compute  the  wavelet  transform  at  scales  smaller  than  1 .  Therefore,  if  the  uniform  Lipschitz 
regularity  is  positive.  Equation  68  implies  that  the  amplitude  of  the  wavelet  transform 
modulus  maxima  should  increase  when  the  scale  increase. 

We  can  define  the  negative  Lipschitz  exponents  for  tempered  distribution  as  [50,  15];  A 
tempered  distribution  is  uniformly  Lipschitz  a  on  ]a,,b[  if  and  only  if  its  primitive  is 
uniformly  Lipschitz  a  +  1  on  ]a,  b[. 

For  example,  the  primitive  of  a  Dirac  function  6{x  —  Xo)  is  a  function  which  is  a  step  edge 
at  Xq.  This  step  edge  is  uniformly  Lipschitz  a  =  0  in  the  neighborhood  xq.  Hence,  a  Dirac 
is  thus  equal  to  —1  in  the  neighborhood  of  Xq  and  its  wavelet  maxima  should  decrease  as 
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the  scale  increase.  This  can  indeed  be  verified  in  Figure  18.  Figure  18  (a)  shows  3  Gaussian 
bumps  with  different  variance  respectively,  1600,20,1.  We  projected  these  images  on  to 
a  second  derivative  of  a  cubic  spline  wavelet  at  scales  2,  4,  8,  16,  32  and  64.  The  result  is 
shown  in  Figure  18  (b).  The  object  at  abscissa  3  produces  wavelet  transform  maxima  at 
decrease  as  the  scale  increases.  Its  distribution  is  locally  equal  to  a  Dirac.  The  wavelet 
transform  maxima  of  the  object  at  abscissa  2  first  increases  then  decreases  as  the  scale 
increases.  We  know  that  as  the  scale  increases,  the  size  of  the  neighborhood  also  increases. 
When  the  neighborhood  increases  to  a  certain  size,  the  regularity  changes.  Therefore,  the 
wavelet  transform  maxima  begin  to  decrease.  Therefore,  the  local  maximum  across  scale  is 
associated  with  the  boundary  of  an  object,  while  the  local  maximum  across  the  shifting 
parameter  indicates  the  location  of  an  object. 

The  relationship  between  local  maximum  across  scale  and  the  boundary  of  an  object  is 
further  demonstrated  in  Figure  19.  The  left-hand  side  of  Figure  19(a)  shows  three  image 
profiles,  containing  synthetic  bumps  of  distinct  width.  We  calculated  the  wavelet  transform 
for  each  image.  At  each  scale,  the  maximum  of  wavelet  coefficients  across  the  shifting 
parameter  was  found.  The  relation  between  scale  and  the  wavelet  maximum  is  shown  on 
the  right-hand  side  of  Figure  19(a).  Note  that  there  is  exactly  one  scale  a*  corresponding 
to  the  maximum  of  wavelet  coefficients  across  both  shifting  and  scale  parameters.  The 
value  a*  is  clearly  a  function  of  the  size  of  the  feature  within  the  cropped  region. 

Scale  Search  for  Tumor  Detection 

In  this  section,  we  define  a  normalized  contrast  and  propose  a  scheme  to  find  the  local 
maximum  of  the  normalized  contrast  across  scale. 

In  the  previous  discussion,  we  argued  that  the  local  maximum  across  scale  is  associated 
with  the  size  of  an  object,  while  the  local  maximum  across  shifting  parameter  indicate  the 
location  of  an  object.  Since  we  are  targeting  tumors  of  a  certain  size,  we  know  a  priori,  the 
expected  range  of  sizes.  We  consider  the  neighborhood  of  a  wavelet  transform  maxima  as  a 
Region  Of  Interest  (ROI).  Since  our  CFAR  wavelet  detector  utilizes  the  difference  between 
the  test  region  and  its  surrounding,  we  would  like  to  maximize  the  difference  between  the 
ROI  and  its  background  per  unit  distance.  For  a  1-D  image,  let  Xi  be  a  wavelet  transform 
maximum.  The  normalized  contrast  is  defined  as 

G  =  —  ^maXjepoo)  ^  j  _ 

The  normalized  contrast  for  a  2-D  image  is  an  average  of  8  directional  maximal  gradients. 
We  would  like  test  this  index  with  a  mathematic  phantom  shown  in  Figure  20  (a).  In 
Figure  21  (b),  we  show  that  the  normalized  contrast  first  increases  then  decreases,  and  the 
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Figure  19:  Best  scale  selection,  (a)  Left:  1-D  pulses  of  different  widths  are  shown.  Right: 
the  maximum  of  wavelet  maxima  across  scale  can  be  used  to  determine  the  width  of  the 
pulse,  (b)  The  results  in  (a)  are  duplicated  with  a  large  amount  of  simulation  noise. 
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Figure  20:  Best  scale  selection:  (a)  The  ideal  phantom,  (b)  wavelet  coefficients  at  the  best 
scale  (34.88),  (c)  wavelet  coefficients  at  scale  24. 

location  of  this  wavelet  maxima  is  at  the  center  of  the  mass.  When  the  scale  becomes  small 
enough,  the  detector  starts  to  focus  on  smaller  structures,  i.e.,  the  boundary  of  the  mass. 
The  wavelet  maximum  moves  toward  the  edge  of  the  mass  as  the  scale  decreases.  And  the 
wavelet  representation  become  a  ring  as  shown  in  20  (c).  The  local  maximum  across  the 
scale  is  the  best  scale  to  detect  that  mass  since  the  difference  of  ROI  to  its  background  is 

maximized. 

Having  proved  the  case  using  an  ideal  phantom,  we  now  apply  our  scale  search  scheme  to 
medium  size  mass  in  the  RMI  phantom  already  shown  in  12  (a).  Figure  22  shows  a  similar 
result  to  the  ideal  phantom  case.  The  noise  in  the  RMI  radiograph  cause  some  small 
fluctuation  in  the  Contrast-Scale  curve.  If  stuck  at  a  small  local  maximum  during  the 
search  for  maximum,  i,  we  use  a  small  force  to  get  away  from  the  small  local  maxima. 
Eventually,  we  will  reach  the  maximum.  We  also  test  our  algorithm  on  some  real 
mammograms,  shown  in  Figure  14,  16  and  17.  The  results  are  shown  in  Figure  23,  24  and 
25.  All  these  searches  find  a  maximum  across  the  scale.  When  we  use  these  best  scale  m 
our  detector,  we  could  use  higher  the  detection  threshold  than  when  we  randomly  pick  a 

scale. 

2.2.5  Summary 

In  this  chapter,  we  derived  a  novel  continuous  scale  wavelet  detector.  The  importance  of 
arbitrary  scale  analysis  was  demonstrated  by  digital  radiographs  of  a  mammography 
phantom  and  digitized  mammograms.  For  the  first  time,  wavelets  and  CFAR  detectors 
were  combined  into  a  single  formulation.  And  we  demonstrated  that  this  powerful  detector 
was  able  to  detect  very  subtle  masses,  which  were  rated  to  be  almost  invisible  by 
radiologist  specializing  in  mammography. 

We  have  put  maxima-across-scale  scheme  [51,  52]  under  the  wavelet  framework.  We  have 
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shown  that  wavelet  maxima  across  scales  can  provide  the  boundary  information  of  an 
object  and  also  indicate  the  location  of  an  object.  By  combining  wavelet  theory  and  the 
geometric  knowledge  of  a  tumor,  we  propose  to  use  a  normalized  contrast  as  an  index  to 
find  the  best  scale.  This  method  is  shown  to  be  effective  in  the  scale  selection  for  tumor 

detection. 
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2.3  Coherence  of  Multiscale  Features  for  Digital  Mammogram 
Enhancement 

2.3.1  Introduction 

Mammograms  can  detect  tumors  that  are  an  eighth  of  an  inch  in  diameter,  while  manual 
examination  usually  fails  to  detect  tumors  smaller  than  a  half-inch.  Reliable  diagnosis  by 
radiographs  of  malignant  breast  disease  depends  on  observing  local  and  distant  changes  in 
tissues  produced  by  disease.  Of  the  visual  clues  of  cancer  sought  by  radiologists,  the 
preliminary  signs  of  masses  and  calcifications  are  most  important  [53,  54].  Unfortunately, 
at  the  early  stages  of  breast  cancer,  these  signs  are  very  subtle  and  varied  in  appearance, 
making  diagnosis  difficult  and  challenging  even  to  specialists  [53,  28]. 

Radiographic  signs  of  cancer  are  related  to  tumor  mass  and  its  density,  size,  shape,  borders 
and  calcification  content.  Extraction  of  these  features  and  enhancement  of  them  can  assist 
radiologists  to  locate  suspicious  areas  more  reliably  [55]. 

Multiscale  representations  based  on  wavelets  have  been  carried  out  for  mammographic 
feature  analysis  [9,  56,  57].  Laine  et  al.  [9]  used  two  overcomplete  multiscale 
representations  for  contrast  enhancement.  Mammograms  were  reconstructed  from 
transform  coefficients  modified  at  each  level  by  nonlinear  weighting  functions.  Qian  et  al. 
[56]  introduced  tree-structured  nonlinear  filters  for  microcalcification  cluster  detection.  An 
image  was  enhanced  by  tree-structured  nonlinear  filters  with  fixed  parameters  and  adaptive 
order  statistic  filters.  Richardson  Jr.  [57]  applied  linear  and  nonlinear  filtering  approaches 
to  the  analysis  of  mammograms.  Here,  a  linear  multiscale  decomposition  was  obtained  via 
a  wavelet  transform;  a  nonlinear  multiscale  decomposition  employed  a  “mean  curvature 
partial  differential  equation”  filter  and  “weighted  majority-minimum  range”  filter.  In 
addition,  Li  et  al.  [58]  extended  a  conventional  multiresolution  wavelet  transform  into  a 
multiresolution  and  multiorientation  wavelet  transform.  They  applied  directional  wavelet 
analysis  to  capture  orientation  information  within  each  mammogram. 

Freeman  and  Adelson  [29]  first  proposed  the  concept  of  steerable  filters  and  applied  it  to 
several  problems  in  the  area  of  computer  vision.  With  a  set  of  “basis  filters”,  one  can 
adaptively  steer  a  filter  along  any  orientation.  Hilbert  transform  pairs  were  constructed  to 
find  a  local  “oriented  energy”  measure  and  dominant  orientation. 

Kass  and  Witkin  [59]  developed  an  algorithm  for  estimating  the  orientation  of  texture 
patterns.  An  orientation  pattern  was  decomposed  into  a  flow  field,  describing  the  direction 
of  anisotropy,  and  a  residual  pattern  obtained  by  describing  the  image  in  a  coordinate 
system  built  from  the  flow  field.  The  texture  orientation  was  estimated  from  Laplacian  of 
Gaussian  filters.  Rao  and  Schunck  [60]  proposed  another  algorithm  based  on  the  gradient 
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of  the  Gaussian.  Their  new  algorithm  incorporated  a  more  sophisticated  scheme  for 
computing  the  coherence  of  the  flow  fleld. 

In  this  chapter,  an  enhancement  algorithm  based  on  overcomplete  multiscale  wavelet 
analysis  is  described.  These  redundant  representations  provide  the  property  of  shift 
invariance  which  is  shown  to  be  advantageous  for  our  analysis.  Features  were  extracted  by 
separable  steerable  filters.  A  coherent  image  and  phase  information  were  then  generated.  A 
nonlinear  function,  integrating  coherent  image  and  phase  information,  was  applied  to  the 
transform  coefficients  at  each  level  of  analysis.  An  enhanced  image  was  obtained  via  an 
inverse  wavelet  transform  of  the  modified  coefficients.  The  novelty  and  advantage  of  this 
algorithm  compared  to  existing  techniques  lies  in  its  detection  of  directional  features  and 
the  lack  of  perturbations  (artifacts)  in  processed  images. 


2.3.2  Background 

In  this  section  we  briefly  describe  the  mathematical  background  and  fundamental  ideas 
used  in  subsequent  sections. 

Wavelet  Transforms 

Wavelet  transforms  have  become  acknowledged  as  useful  tools  for  many  applications  in 
signal  processing.  A  function  i;{x)  is  said  to  be  a  wavelet  if  and  only  if  its  Fourier 
transform  '0(^)  satisfies  the  admissibility  condition 


This  condition  implies  that 


(70) 


This  means  that  ■0(x)  will  have  at  least  some  oscillations. 

Wavelets  constitute  a  family  of  functions  derived  from  one  single  function  {mother 
wavelet)  by  dilations  and  translations 
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where  a  E  R'^,b  e  R.  The  idea  of  the  wavelet  transform  is  to  represent  any  function  f{t)  as 
a  superposition  of  wavelets.  The  continuous  wavelet  transform  is  defined  as 

Wf{a,b)  =  {f,A,b)  =  -JqI 
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where  ip  denotes  the  complex  conjugate  of  A  function  can  be  reconstructed  from  its 
wavelet  transform  by  means  of  the  “resolution  of  identity”  formula 
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In  the  application  of  digital  mammography,  one  prefers  to  express  /  as  a  discrete 
superposition.  Let  us  discretize  the  translation  and  dilation  parameters  of  the  wavelet  in 
Equation  (71): 


_  m 

=  “o  ^  -nbo), 


(74) 


where  a  =  a^, 6  =  nboa^,  with  m,ne  Z,  and  oq  >  1, 6o  /  0.  On  this  discrete  grid,  the 
wavelet  transform  is  simply 

-m.  r+°°  ,  , 

IE/(m, n)  =  ao^  /  f{x)ip{aQ'”‘x  -  nbci)dx.  (75) 

J  — OO 

The  original  signal  can  be  approximated  as  linear  combinations  of  the  wavelet  bases, 

f{x)  ^^Wf{m,n)ipmA^)-  (76) 

m^n 

One  popular  discretization  is  to  choose  Oq  =  2,  6o  =  1, 

-  n),  (77) 


which  are  called  dyadic  wavelets. 

Steerable  Filters 

A  function  f{x,y)  is  called  “steerable”  if  it  can  be  expressed  as  a  linear  combination  of 
rotated  versions  of  itself.  The  fundamental  idea  of  steerable  filters  is  to  apply  “basis  filters” 
which  correspond  to  a  fixed  set  of  orientations  and  interpolate  between  each  discrete 
response.  Thus,  one  must  decide  the  number  of  “basis  filters”  and  the  corresponding 
interpolation  functions.  Let  9i  be  the  angle  of  ith  basis  filter  and  ki{d)  denote  an 
interpolation  function.  As  defined  in  [29]  a  steering  constraint  may  be  formulated  by 

M 

=  p8) 

where  M  is  the  number  of  basis  functions  required  to  steer  a  function  f^^{x,  y). 

Hereafter,  it  will  be  more  convenient  to  work  in  polar  coordinates  r  =  ^x"^  +  y"^  and 
(p  =  arg(a;,  y).  Let  /  be  any  function  that  can  be  expressed  as  a  Fourier  series  in  polar 
angle  (p: 

L 

/(r,0)=  ^  (79) 

n——L 
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where  j  =  and  L  is  the  length  of  the  coefficients. 

Local  Energy 

Based  on  physiological  evidence,  a  local  energy  operator  was  first  proposed  by  Morrone  and 
Owens  [61].  Rather  than  considering  differential  properties,  they  studied  the  properties  of 
the  Fourier  expansion  of  an  intensity  function.  They  demonstrated  that  the  maxima  of  an 
energy  function  are  coincident  with  the  points  of  maximum  phase  congruency  and  visual 
features. 

Let  IG{0)  and  IG^{9)  be  the  outputs  of  constituent  filter  G{e)  and  its  quadrature  pair 
G^{6),  respectively.  The  oriented  energy  function  can  be  expressed  by  the  squared  outputs 
of  a  quadrature  pair  of  filters  steered  to  an  angle  9,  that  is, 

E{9)  =  [IG{9)f  +  [IG^{9)f.  (80) 

Since  IG{9)  and  IG^{9)  can  be  expressed  as  a  sum  of  basis  filter  outputs  times 
interpolation  functions,  following  [29],  Equation  (80)  can  be  simplified  to  a  Fourier  series  in 
angular  form  using  substitutions  obtained  by  the  trigonometric  identity: 

E{9)  =  ao  +  Oi  cos(20)  +  a2  sin(20)  +  [  higher  order  terms •  •  •  ]. 

The  lowest  frequency  terms  of  this  identity  were  used  to  approximate  the  dominant 
directions  A  and  their  associated  magnitudes  M, 

A  =  arctan[oi,  a2]/2,  (81) 

M  =  ^al  +  aj.  (82) 

Figure  26  shows  the  oriented  energy  map  of  a  two-dimensional  Gaussian  function 

f{x,y)=e-^^^^y^K 

The  test  function  is  shown  in  Figure  26(a).  Using  a  second  derivative  of  Gaussian  function 
and  its  Hilbert  transform  to  measure  oriented  energy,  we  plot  the  direction  and  magnitude 
of  each  position  on  a  map  where  the  arrow  indicates  the  direction  and  the  length  is  the 
magnitude,  which  is  shown  in  Figure  26(b).  The  edges  were  detected  at  the  inflection 
points  of  the  surface.  The  magnitudes  of  the  energy  function  are  maxima,  therefore  the 
lengths  are  longer  than  those  of  other  points  of  the  function. 

Note  that  Figure  26(b)  does  not  show  the  correct  phases  in  the  upper  half  of  the  map.  All 
the  arrows  are  pointing  outward  instead  of  inward.  There  is  a  180°  phase  shift.  This  is 
because  the  oriented  energy  is  obtained  from  the  squared  outputs  of  quadrature  pairs.  We 
applied  a  simple  algorithm  to  correct  this  phase  shift.  First,  we  divide  the  space  into  8 
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(a)  A  test  function. 


(b)  The  orientation  and  magnitude  map  of  a  test  function. 


Figure  26:  An  example  of  oriented  energy  map  construction 


directions  and  approximate  the  phase  at  each  point  to  exactly  one  of  these  discrete 
orientations.  Next,  the  neighbors  of  each  point  along  the  approximated  direction  are 
examined.  If  the  intensities  along  that  direction  are  increasing,  the  phase  is  not  changed; 
otherwise,  we  correct  the  phase  by  a  180°  shift.  A  sample  corrected  oriented  energy  map  is 
shown  in  Figure  27. 

Measure  of  Coherence 

Texture  plays  an  important  role  in  many  machine  vision  and  image  processing  tasks 
including  surface  inspection,  scene  classification,  surface  orientation  and  shape 
determination  [62].  Texture  patterns  may  be  characterized  by  extracting  measurements 
that  quantify  the  nature  and  directions  of  patterns.  Most  breast  carcinomas  have  the 
appearance  of  stellate  lesions  consisting  of  a  central  mass  surrounded  by  radiating  spicules 
[53].  The  spicules  radiate  outward  in  all  directions  and  vary  in  length.  This  provides  an 
important  cue  for  the  early  detection  of  such  lesions. 

Much  attention  has  been  given  to  the  notion  of  decomposing  an  intensity  image  into 
intrinsic  images  to  extract  meaningful  information  [63,  64].  Such  intrinsic  properties 
represent  basic  components  of  the  image  formation  process  and  therefore  can  reveal 
features  “hidden”  inside  an  image.  The  information  they  provide  is  beyond  the  intensity 
image  alone. 

Rao  and  Schunck  [60]  defined  the  orientation  field  of  a  texture  image  to  consist  of  two 
images  —  an  angle  image  and  a  coherence  image.  The  angle  image  denotes  the  dominant 
local  orientation  at  each  point  while  the  coherence  image  represents  the  degree  of 
anisotropy  at  each  point.  They  strongly  advocated  the  use  of  angle  and  coherence  images 
as  intrinsic  images.  In  this  paper,  we  investigated  the  efficiency  of  these  two 
representations  to  capture  and  enhance  features  of  importance  to  mammography. 

2.3.3  Methodology 

Our  algorithm  consists  of  four  steps:  (1)  construction  of  an  overcomplete  multiscale 
representation,  (2)  processing  via  separable  steerable  filter  analysis,  (3)  computing 
coherence  measures  and  (4)  application  of  nonlinear  operators,  as  shown  in  Figure  28  [65]. 
We  next  describe  each  block  in  the  sections  below. 

(1)  Overcomplete  Multiscale  Representations:  Wavelet  transforms,  owing  to  their 
localization  characteristics,  are  powerful  tools  of  analysis  for  many  signal  and  image 
processing  applications.  Multiscale  analysis  can  extract  features  at  distinct  scales  and 
provide  local  information  often  hidden  in  an  original  mammogram.  One  major  drawback  of 
wavelet  transforms  is  their  lack  of  translation  invariance  [45],  making  the  content  of 
wavelet  subbands  unstable  under  the  translations  of  an  input  signal.  In  our  algorithm,  a 
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Figure  28:  Block  diagram  of  the  proposed  algorithm. 

digitized  mammogram  was  decomposed  using  a  fast  wavelet  transform  algorithm  (FWT) 
[15].  In  order  to  obtain  wavelet  coefficients  at  each  level  without  downsampling,  an 
undecimated  ^^algorithme  a  trous"  (algorithm  with  holes)  [21,  66]  was  implemented.  In  the 
spatial  domain,  the  redundant  representation  corresponds  to  no  aliasing. 

Let  denote  an  original  image  and  D*/  be  obtained  by  inserting  2*  —  1  zeros  between 
every  pair  of  the  coefficients  representing  /.  (!>*/)»  and  {D^f)y  stand  for  carrying  out 
convolution  operations  with  the  filter  D^f  along  x  and  y  directions,  respectively.  The 
decomposition  and  reconstruction  equations  at  level  i  are  as  follows: 

Decomposition: 


Si+i 

=  *{D^h)^*{D%)y, 

=  S^*{D^g)x. 

YY  y 

=  SU{D^g)y. 

(83) 

Reconstruction: 

*  (D^k),  *  {DH)y  +  *  (DH),  *  {D^k)y  +  5*+'  *  (D^h)^  *  {D%y,  (84) 

where  indicates  discrete  convolution,  and  h,g,k,  and  I  are  filters  whose  Fourier 
transforms  {H{oj),G{u)),  K{uj),  and  L{u),  respectively)  satisfy  [15] 


G{to)K{u)  +  \H{oj)\^  =  l, 
l  +  \H{u)\^ 

L{^)  = - ^ - • 


(85) 
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Figure  29:  Block  diagram  of  a  fast  wavelet  transform  at  level  i.  (a)  Decomposition,  (b) 
Reconstruction. 

The  block  diagrams  of  decomposition  and  reconstruction  operations  at  level  i  were  shown 
in  Figure  29(a)  and  (b),  respectively. 

(2)  Separable  Steerable  Filter  Analysis:  A  filter  is  called  “steerable”  if  the  filter  at  an 
arbitrary  orientation  can  be  expressed  as  a  linear  combination  of  a  set  of  basis  filters, 
generated  from  rotations  of  a  single  kernel  [29].  Steerable  filters  [29,  67,  68],  which  can  be 
adaptively  adjusted  to  arbitrary  orientation,  were  used  to  detect  stellate  patterns  of 
spicules  and  locate  feature  orientations  more  precisely.  As  pointed  out  by  [69],  the 
separability  property  of  the  filters  can  speed  up  computations  considerably  when  convolved 
with  a  large  image  matrix.  In  our  algorithm,  we  used  three  basis  functions  as  steerable 
filters.  The  x-y  separable  steerable  approximations  of  filter  kernels  were  generated  by 
Singular  Value  Decomposition  (SVD)  [29,  69].  Using  a  set  of  separable  steerable  filters,  the 
magnitude  (M*)  and  associated  dominant  directions  (A*)  of  local  energy  were  determined 
by  the  basis  functions  and  their  Hilbert  transforms  [29,  70,  61,  71].  Figure  30  illustrates  the 
detailed  procedure  of  steerable  filter  processing.  We  applied  separable  steerable  filters  to  an 
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original  image  5**^  and  a  low-pass  filtered  image  at  each  level  S^,i  =  1 , . . .  n,  to  capture 
salient  multiscale  features.  Discrete  realizations  of  the  three  basis  functions  were  applied  in 
our  implementation  algorithm.  An  approximation  of  Hilbert  transform  of  basis  function 
was  thereafter  obtained. 

(3)  Coherence  Measures:  A  coherence  map  is  an  image  showing  a  local  measure  of  the 
degree  of  anisotropy  of  flow  [59,  60].  If  the  orientations  of  a  texture  pattern  at  any  point 
{xi,yi)  are  coherent,  then  magnitude  and  phase  information  are  important  and  should  be 
emphasized.  Conversely,  if  the  orientations  are  not  coherent,  the  magnitude  and  phase 
information  can  be  neglected  or  attenuated.  Kass  and  Witkin  [59]  suggested  a  simple  way 
of  measuring  strength  of  coherence  by  finding  the  ratio 


. .  _  \W{j,k)J{j,k)\ 

{w{j,k)\j{j,k)\y 


(86) 


where  W{j,k)  denotes  a  local  weighting  function  with  unit  integral,  J{j,k)  denotes  the 
squared  gradient  vector  at  {j,k),  and  |  •  |  denotes  absolute  value. 

An  alternative  measure  of  coherence  was  proposed  by  Rao  and  Schunck  [60]  and  was 
obtained  by  weighting  the  energy  with  the  normalized  projection  of  energy  within  a 
specified  window  (W)  onto  the  central  point  (j,  k)  of  a  window.  The  coherence  measure 
(C‘)  was  expressed  as 


C^{j,k)=Myj,k) 


^2  \Mym,n)cos{Ayj,k) 

(m,n)GW 

^2  M^{m,n) 

(m,7i)6W 


A*(m,  n))| 
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(87) 


where  M^j^k)  and  Ayj,k)  denote  energy  and  phase  of  point  {j,k)  at  level  i,  respectively. 
This  coherence  measure  incorporates  the  gradient  magnitude  and  hence  places  more  weight 
on  regions  that  have  higher  visual  contrast.  Let  us  examine  more  closely  how  the  coherence 
measure  works.  Suppose  that  there  is  a  region  somewhere  in  an  image,  as  shown  in 
Figure  31.  Pixel  A  is  on  a  vertical  edge  line  because  its  magnitude  is  maximum  along  the 
direction  of  its  phase.  (Actually,  there  is  an  edge  along  the  vertical  direction  of  pixel  A.) 
The  coherence  at  pixel  A  is  measured  inside  Window  1  which  surrounds  pixel  A.  Since  the 
pixels  inside  Window  1  have  the  same  phases  as  pixel  A,  the  coherence  of  pixel  A  is  equal 
to  the  magnitude  at  this  pixel.  While  the  pixels,  surrounding  pixel  B,  inside  Window  2 
point  to  distinct  directions,  the  projections  of  these  pixels  on  pixel  B  will  be  balanced.  The 
coherence  at  pixel  B  is  reduced.  This  pixel  can  be  considered  a  disturbance. 
Experimentally,  this  method  out-performed  previous  measures  applied  to  similar  data  [60]. 
In  our  algorithm,  we  implemented  this  approach  to  measure  each  coherence  map.  In  order 
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Figure  30:  Diagram  of  steerable  filter  processing 


Figure  31:  An  example  showing  the  effectiveness  of  coherence  measure. 

to  capture  features  at  each  level,  we  changed  the  window  size  at  each  distinct  level.  At  the 
finest  level,  a  5  x  5  window  was  used.  While  for  the  coarser  levels,  larger  windows  (e.g. 

7  X  7  or  9  X  9)  were  used.  For  simplicity  of  implementation,  an  odd  window  size  was 
selected  at  each  level.  Figure  32  shows  an  example  of  a  5  x  5  window  used  to  carry  out  this 
operation.  The  measure  of  coherence  C*  was  obtained  from  Equation  (87).  The 
combination  of  coherence  and  orientation  structure  was  better  able  to  extract  the  salient 
features  of  spiculated  lesions.  In  overview,  a  schematic  diagram  of  processing  Steps  1-3  is 
shown  in  Figure  33. 

(4)  Nonlinear  Operators:  We  have  now  computed  all  the  information  needed  complete 
algorithm.  A  nonlinear  operation  is  next  applied  within  each  level  to  precisely  modify 
transform  coefficients.  This  operation  integrates  both  coherence  map  and  phase 
information,  as  described  below: 

A.  Modification  from  Coherence  Maps: 

Let  C^{j,  k)  denote  the  coherence  measure  of  point  (j,  k)  at  some  level  i.  First,  we  find  the 
local  maxima  of  the  coherence  map  at  each  level.  These  maxima  correspond  to  features  in 
an  image  and  their  propagations  at  distinct  levels  of  analysis.  Therefore  they  are  likely 
locations  to  “boost”  or  amplify.  We  call  these  maps  (Cj)  feature  maps.  Then  the  maxima 
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Figure  32:  Illustration  of  measuring  the  coherence  of  an  image. 


or  square  roots  of  these  maxima  are  mapped  into  256  scales,  depending  on  the  threshold 
value  of  the  scale  (5^*)  at  the  first  level  (We  will  describe  how  to  obtain  this  value  later). 
That  is, 


255^-i — Lsan 


^  fyrnin 


255 


if  Si  <  St, 
otherwise, 


where  and  are  maximum  and  minimum  of  C},  respectively.  St  is  a 

prespecified  threshold  value.  We  compute  the  square  root  of  the  maxima  because  the 
difference  between  the  maximal  values  might  be  large  (especially  at  the  first  two  levels). 
This  operation  prevents  most  values  mapping  to  a  small  scale. 

We  then  construct  a  histogram  of  square  roots  of  local  maxima  and  accumulate  the  number 
of  mappings  from  1  up  to  some  scale  {Sj),  so  that  the  accumulated  number  is  over  99%  of 
the  total  local  maxima.  The  corresponding  value  of  that  scale  is  then  chosen  as  the 
threshold  value  (C'),J  for  that  level. 

Modifications  of  coherence  maps  were  obtained  by  a  nonlinear  function  expressed  as 


CitU,k) 


if  C\j,k)<C%{j,k), 
otherwise, 


(88) 


where  Cl^eanU^  k)  is  the  mean  value  of  the  figure  map  at  each  z-th  level.  For  an  image  with 
a  black  background,  the  background  pixels  are  not  be  counted  in  computing  the  mean 
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Figure  33;  Overview  of  processing  for  steps  1-3. 
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value.  The  square  root  values  were  computed  so  that  we  would  not  over-emphasize  large 
local  maxima  and  ignore  the  small  local  maxima. 

B.  Modification  from  Phase: 

Phase  information  is  important  to  distinctly  characterize  oriented  texture.  Therefore  we 
did  not  neglect  its  contribution  in  the  modification  of  coefficients.  We  applied  a  sinusoidal 
weighting  to  phase  information.  The  detail  sub-bands  of  wavelet  coefficients  obtained  in 
Step  1  included  two  components:  the  component  along  the  x  direction  and  the  component 
along  the  y  direction.  The  x  component  was  obtained  by  high-pass  filtering  along  the  x 
direction,  hence  mostly  vertical  features  within  the  mammogram  were  detected.  We 
emphasized  the  points  whose  dominant  orientations  were  near  0  and  tt  (with  respect  to  the 
vertical  axis).  Thus,  the  modification  from  phase  information  was 

KiodU^  k)  =  I  cos{A\j,  k))\.  (89) 

The  y  component  was  obtained  by  high-pass  filtering  along  the  y  direction,  detecting 
horizontal  features.  We  emphasized  the  points  whose  dominant  orientations  were  near  | 
and  ^  (with  respect  to  the  vertical  axis).  Transform  coefficients  were  modified  by  phase 
information  by 


^modih  k)  =  I  sin(^*(j,  k))\.  (90) 

The  modifications  from  coherence  map  and  phase  at  level  i  —  1  were  conibined  to  adjust 
the  wavelet  coefficients  at  level  i.  The  resultant  modification  was  therefore 


rpi-l  _ 
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where  T*  was  a  constant  at  each  level.  A  schematic  diagram  of  this  nonlinear  operator  at 
level  i  is  shown  in  Figure  34.  Finally,  the  modified  coefficients  were  then  used  to 
reconstruct,  via  an  inverse  fast  wavelet  transform,  an  enhanced  visualization  of 
mammographic  features. 


2.3.4  Experimental  Results  and  Discussion 

In  this  section,  we  present  some  examples  of  processed  results.  First,  we  show  a 
one-dimensional  example  which  does  not  take  advantages  of  oriented  information.  We 
perform  an  overcomplete  wavelet  transform  to  the  signal  and  decompose  it  into  five  levels. 
The  enhanced  signal  and  the  original  signal  are  plotted  together  in  Figure  35  to  show  the 
benefit  of  our  algorithm.  The  transition  regions  in  the  original  signal  are  clearly  enhanced. 
Three  examples  of  malignant  lesions  with  distinct  radiographic  signs  of  cancer  were 
processed.  The  images  were  of  matrix  size  512  x  512.  For  each  case,  both  global  and 
regions  of  interest  (ROI)  are  shown  along  with  the  corresponding  enhanced  ROI  image. 
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Microcalcifications 

Figure  36  shows  a  sample  mammogram  (mamOOSrml)  with  microcalcification  clusters.  The 
original  mammogram  is  shown  in  Figure  36(a).  Figure  36(b)  shows  an  original  suspicious 
area  of  the  mammogram.  After  enhancement,  clusters  of  calcifications  clearly  appear  in  the 
center  of  the  image. 

Stellate  lesions 

A  mammogram  (mam041rml)  with  a  stellate  lesion  is  shown  in  Figure  37.  Figure  37(a) 
shows  the  original  image.  Figure  37(b)  presents  an  original  digital  radiograph  with  a 
partially  obscured  irregular  mass  in  the  center  of  the  image  matrix.  After  applying  our 
algorithm  to  this  image,  the  enhanced  image  makes  obvious  spiculated  lesions  around  the 
mass.  The  mass  itself  was  also  enhanced,  as  shown  in  Figure  37(c).  These  clear  spiculated 
patterns  suggested  radiologists  that  this  mass  is  more  likely  malignant. 

Masses 

A  mammogram  (mam0041cc)  with  a  mass  tumor  is  shown  in  Figure  38.  The  craniocaudal 
view  of  the  left  breast  shown  in  Figure  38(b)  contains  an  irregular  spiculated  mass  in 
retroglandular  fat.  The  enhanced  ROC  shown  in  Figure  38(c)  better  delineates  the  margins 
of  the  mass. 

Figures  39(a)  and  (b)  show  five  levels  of  wavelet  coefficients  along  the  x  and  y  directions, 
respectively.  The  corresponding  coherence  maps  are  shown  in  Figure  40(a).  Feature  maps 
are  shown  in  Figure  40(b).  The  histogram  of  square  roots  of  local  maxima  are  constructed 
in  Figure  40(c)  and  threshold  scales  were  determined  as  24,  32,  81,  109,  and  123  for  five 
levels.  The  wavelet  coefficients  were  then  modified  prior  to  reconstruction. 

To  validate  our  enhancement  techniques,  three  mathematical  models  of  phantoms  were 
constructed.  The  models  included  features  of  interest  in  mammographic  imaging,  such  as 
masses,  microcalcifications,  and  spicular  lesions,  as  shown  in  Figure  41.  These  phantoms 
were  blended  into  a  normal  mammogram  to  form  an  experimented  database.  Note  that  the 
positions  of  these  mathematical  models  are  different  in  each  image.  We  constructed  twenty 
cases  in  our  testing  database.  The  images  were  of  matrix  size  1024  x  1024. 

Many  techniques  of  contrast  manipulation  and  modification  have  been  developed  within 
the  field  of  digital  image  processing.  However,  the  measurement  and  evaluation  of  contrast 
and  contrast  changes  in  arbitrary  images  are  not  uniquely  defined  in  the  literature  [72,  73]. 
In  this  paper,  we  adopt  a  definition  introduced  by  Morrow  et  al.  [10] 


a  = 


Bo-B 
Bo  +  B’ 


(92) 


where  B  is  the  mean  gray-level  value  of  a  particular  object  in  the  image  (foreground),  and 
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Figure  41:  Mathematic  phantoms. 


Bq  is  the  mean  gray-level  value  of  a  surrounding  region  (background). 

Based  on  this  definition,  we  applied  the  proposed  algorithm  to  each  blended  image  and 
computed  distinct  contrast  values.  A  quantitative  measure  of  contrast  improvement  was 
defined  by  a  contrast  improvement  index  (CII)  [9,  74]: 


CII  = 


^processed 
^original  ’ 


where  Qonginai  contrast  values  for  a  region  of  interest  in  the  processed 

and  original  images,  respectively. 

We  compared  our  results  with  those  of  three  existing  contrast  enhancement  techniques: 
contrast  stretching  (CS),  unsharp  masking  (UM),  and  histogram  equalization  (HE),  as 
shown  in  Table  8-10.  From  these  results,  we  observe  that  the  values  of  CII  depend  on  the 
image  itself,  the  types  of  features  and  the  surroundings  context  of  the  features.  The 
contrast  stretching  and  unsharp  masking  methods  consistently  had  a  little  improvement 
over  the  original  image,  while  the  histogram  equalization  approach  had  varied  performance. 
Our  algorithm  consistently  outperformed  these  other  methods.  Note  that  since  masses 
have  larger  regions  than  those  of  microcalcifications  and  spicular  lesions.  Since  we 
emphasized  features  at  edges  the  mean  values  of  the  masses  were  less  improved. 


Table  8:  Contrast  Improvement  Index  for  masses  obtained  by  Contrast  Stretching  (CS), 
Unsharp  Masking  (UM),  Histogram  Equalization  (HE),  and  proposed  Coherent  Features 

(CF)  methods. 


Case 

CS 

UM 

HE 

CF 

1 

1.2321 

1.1538 

0.3934 

3.4121 

2 

1.2720 

1.1627 

0.5430 

3.5911 

3 

1.2453 

1.1370 

1.1606 

4.1997 

4 

1.2218 

1.1109 

0.6830 

3.5410 

5 

1.2313 

1.1125 

0.2755 

3.7025 

6 

1.2590 

1.2214 

2.3413 

4.2358 

7 

1.2458 

1.2080 

3.0524 

4.9465 

8 

1.2374 

1.0895 

0.9045 

4.2415 

9 

1.3271 

1.1393 

1.7784 

6.1136 

10 

1.2236 

1.1569 

0.8641 

3.9763 

11 

1.2622 

1.0887 

1.2967 

3.6320 

12 

1.2485 

1.1150 

3.7718 

5.9985 

13 

1.2543 

1.1119 

1.1636 

3.1344 

14 

1.2742 

1.1809 

2.5198 

4.7512 

15 

1.2332 

1.1084 

2.0114 

3.8008 

16 

1.2316 

1.1533 

0.9142 

4.1431 

17 

1.2071 

1.1520 

0.1207 

3.0553 

18 

1.2233 

1.1074 

2.4443 

3.9690 

19 

1.2275 

1.1273 

2.0329 

3.1121 

20 

1.2614 

1.1678 

1.9330 

3.8542 

2.3.5  Summary 

An  enhancement  algorithm  relying  on  multiscale  wavelet  analysis  and  extracting  oriented 
information  at  each  scale  of  analysis  was  investigated.  The  evolution  of  wavelet  coefficients 
across  scales  characterized  the  local  shape  of  irregular  structures.  Using  oriented 
information  to  detect  the  features  proved  to  be  an  effective  method  of  enhancing  complex 
and  subtle  structures  of  the  breast.  Steerable  filters,  rotated  at  arbitrary  orientations 
reliably  found  visual  cues  within  each  spatial-frequency  sub-band  of  an  image.  Coherence 
measure  and  dominant  orientation  clearly  helped  to  discriminate  features  from  complex 
surrounding  tissue  in  mammograms. 


95 


Table  9;  Contrast  Improvement  Index  for  microcalcifications  obtained  by  Contrast  Stretch¬ 
ing  (CS),  Unsharp  Masking  (UM),  Histogram  Equalization  (HE),  and  proposed  Coherent 
Features  (CF)  methods. 


Case 

CS 

UM 

HE 

CF 

1 

1.2510 

2.0225 

3.0738 

8.0812 

2 

1.2855 

1.9996 

1.5135 

8.8415 

3 

1.2811 

1.9945 

2.1976 

9.0558 

4 

1.2204 

1.9953 

0.6681 

6.1101 

5 

1.2322 

1.9159 

0.3781 

6.5205 

6 

1.2703 

1.9731 

1.3389 

7.8585 

7 

1.2335 

1.9508 

0.6793 

6.5333 

8 

1.2394 

2.0097 

1.1124 

7.7731 

9 

1.2636 

1.7277 

1.7397 

6.5394 

10 

1.2280 

2.1050 

0.8854 

5.8531 

11 

1.3195 

1.7763 

1.5760 

6.1981 

12 

1.2354 

1.9701 

2.0045 

6.8729 

13 

1.2938 

2.1740 

3.4047 

6.1026 

14 

1.2404 

1.9956 

0.1757 

7.3612 

15 

1.2840 

2.1781 

1.0409 

5.8868 

16 

1.2452 

1.9345 

2.1372 

6.6954 

17 

1.2438 

1.9107 

0.9802 

6.4960 

18 

1.2246 

1.9772 

2.3434 

6.5619 

19 

1.2299 

1.9302 

2.1773 

6.3651 

20 

1.2308 

1.9459 

1.0755 

5.9010 
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Table  10:  Contrast  Improvement  Index  for  spicular  lesions  obtained  by  Contrast  Stretch¬ 
ing  (CS),  Unsharp  Masking  (UM),  Histogram  Equalization  (HE),  and  proposed  Coherent 
Features  (CF)  methods. 


Case 

CS 

UM 

HE 

CF 

1 

1.2576 

1.5373 

3.5718 

6.4958 

2 

1.3016 

1.6564 

1.5741 

8.9896 

3 

1.2709 

1.6748 

1.5917 

9.0420 

4 

1.2647 

2.1918 

2.6259 

14.8793 

5 

1.2397 

1.5258 

1.0526 

6.4616 

6 

1.2508 

1.5366 

1.5741 

6.4331 

7 

1.2432 

1.5318 

2.1020 

6.4341 

8 

1.2837 

1.7600 

1.1114 

11.1680 

9 

1.3500 

1.6194 

1.0969 

10.5974 

10 

1.2203 

1.6627 

0.6285 

7.1692 

11 

1.2498 

1.6466 

1.7991 

7.1298 

12 

1.2654 

1.6334 

2.1535 

8.9059 

13 

1.2835 

1.6544 

1.0977 

5.4763 

14 

1.2495 

1.6130 

0.6353 

7.8101 

15 

1.2860 

1.7091 

1.0337 

5.3742 

16 

1.2332 

1.6502 

1.1970 

7.2288 

17 

1.3066 

1.4405 

1.6293 

6.3881 

18 

1.2135 

1.5815 

0.8470 

6.2905 

19 

1.2527 

1.5995 

3.1518 

7.1808 

20 

1.2478 

1.6107 

3.3187 

6.1226 
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2.4  A  Lifting  Algorithm  for  Overcomplete  Wavelet  Representations  for 
Interactive  Feature  Analysis 

2.4.1  Introduction 

In  this  section  we  first  review  the  basic  idea  behind  the  lifting  scheme  introduced  by  Wim 
Sweldens  [75]  [76]  as  a  flexible  tool  for  constructing  compactly  supported  biorthogonal  wavelets. 
We  study  in  Section  2.4.2  the  properties  of  second  generation  wavelets  that  emerged  from  the 
need  of  more  general  basis  functions  in  some  applications  on  general  domains.  In  Section  2.4.3, 
we  review  the  biorthogonal  wavelets  and  duals  to  see  how  they  are  related  each  other  to  satisfy  a 
perfect  reconstmction  condition  in  a  filter  bank  coding  system.  In  Section  2.4.4  we  describe  the 
mathematical  background  of  the  lifting  scheme  and  its  definition.  We  then  discuss  in  Section  2.4.5 
the  relation  between  interpolating  scaling  functions  and  the  lifting  scheme.  In  Section  2.4.6  we 
show  how  this  lifting  scheme  can  be  generalized  to  construct  any  wavelet  or  wavelet  transforms 
with  a  finite  number  of  lifting  steps  [77].  Finally,  in  Section  2.4.7  we  introduce  an  overcomplete 
lifting  scheme  that  can  enhance  the  performance  of  the  wavelet  processing  via  overcomplete 
representations. 

2.4.2  Second  Generation  Wavelets 

Wavelets  ^  are  traditionally  defined  by  translates  and  dilates  of  one  particular  function, 
a  mother  wavelet  v|/ ,  such  that  jfx)  =  \^{a’x-k) .  These  wavelets  are  referred  to  as  first 
generation  wavelets.  The  properties  of  first  generation  wavelets  are  the  following:  They  form  a 
Riesz  basis  for  L2(R) ,  so  that  a  function /in  L2(R)  can  be  represented  with  the  wavelet  basis 
as 

/  ==  Ijt.k'fj.K 

j,k 

where  the  coefficients  yjj.  are  computed  by  =  </,  vji  is  a  dual  wavelet  of  \(/.  The 
wavelets  and  their  duals  are  well  localized  in  space  and  frequency.  In  other  words,  the  wavelets 
and  duals  have  compact  supports,  which  are  smooth,  i.e.  they  decay  towards  high  frequencies,  and 
have  vanishing  polynomial  moments,  which  decay  toward  low  frequency.  Wavelets  are  also 
utilized  in  a  framework  of  multiresolution  analysis.  Finally,  they  mostly  rely  on  the  Fourier 
transform  as  a  basic  instrumental  tool  since  the  operations  of  the  wavelet  family  defined  by 
translates  and  dilates  become  algebraic  operations  in  Fourier  space. 

On  top  of  the  properties  of  the  traditional  wavelets,  second  generation  wavelets  have 
emerged  from  the  need  of  more  generalized  analyzing  functions  which  are  not  necessarily 
translates  and  dilates  of  one  basis  function.  Some  applications  on  general  domains  require 
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wavelets  that  are  defined  on  arbitrary,  possibly  non-smooth,  domains  of  R".  A  special  case  is  the 
construction  of  wavelets  on  bounded  domains  without  introducing  end  effects  near  boundaries. 

We  also  need  wavelets  to  analyze  data  on  curves,  surfaces  or  manifolds  independent  of 
parametrization.  There  are  also  problems  in  real  life  to  analyze  irregularly  sampled  data,  which 
can  not  be  achieved  by  first  generation  wavelets  and  traditional  wavelet  transform  algorithms. 

Because  of  the  non-translation  invariant  property  of  its  basis  fimctions,  second  generation 
wavelets  can  no  longer  use  the  Fourier  transform  as  its  construction  tool.  A  question  arises:  How 
may  second  generation  wavelets  be  constructed?  The  lifting  scheme  introduced  by  Wim  Sweldens 
[75]  [76]  provides  a  profound  answer.  The  scheme  was  initially  introduced  as  a  tool  in 
constructing  biorthogonal  wavelets.  Daubechies  and  Sweldens  [77]  later  generalized  the  lifting 
scheme  for  the  construction  of  any  wavelet  by  factoring  wavelet  transforms  through  finite  lifting 
steps.  In  the  next  section  we  first  study  the  properties  of  the  biorthogonal  wavelets  and  duals 
which  is  the  first  footstone  in  developing  the  lifting  scheme. 

2.4.3  Biorthogonal  Wavelets  and  Duals 

The  orthogonality  property  in  wavelets  imposes  a  strong  limitation  in  the  construction  of 
wavelets.  The  Haar  wavelet  is  the  only  real-valued  wavelet  which  is  compactly  supported, 
orthogonal  and  symmetric  [78].  The  scaling  ftmction  (p  e  L2 ,  and  the  wavelet  function  v|/  e  L2 
both  contribute  to  a  multiscale  analysis  and  must  satisfy  the  refinement  relations  in  the  sense  that, 

(p(x)  =  2^h^.(^(2x—k') , 

k  (2.1) 

x|/(x)  =  . 

k 

The  set  of  translated  scale  functions  {(p(x-^r)|A:  e  Z}  form  a  Riesz  basis.  The  refinement 
relations  are  also  applicable  to  their  duals:  9  and  xj) ,  with  coefficients  h  and  respectively, 
such  that 


9(x)  =  —  , 

k  (2.2) 

V(^)  =  22]g^9(2x-A:)  . 

k 

These  duals  also  fit  into  a  framework  of  multiresolution  analysis.  The  duals  are  biorthogonal  to 
their  primals  in  the  sense  that. 
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(2.3) 


<9(x),  v|i(x-A:))  =  {^(x),  (p(x-k))  =  0  , 

<9(x),(p(jc-^)>  =  <vi/(x),  \|/(x-A:))  = 

The  biorthogonality  condition  (2.3)  is  equivalent  to 

^9(C0  +  2^7t)(p(G)  +  2A:7t)  =  ^\j/(co  +  2^7i)y(co  +  2A:71)  =  1, 

*  ^  ^  ^  _  (2.4) 

^\j/(a)  +  2A:7i)(p(G)  +  2A:Tc)  =  ^(p((»  +  2A:7:)\j/(c()  +  2^7c)  =  0. 

k  k 


Combining  the  condition  (2.3)  or  (2.4)  and  the  refinement  relation  (2.1),  we  can  derive  a 
biorthogonal  condition  for  FIR  filters,  h,g,h,  and  g ,  as 

h(G))h((s))  +  h(Gy  +  7t)h(G)  +  %)  =  1 
g(a))g(co)  +  g(co  +  7i)g(co  +  7i)  =  1 
g((0)A(a))+g(CO  +  7l)/!(CO  +  7l)  =  0 
^(co)g(a))  +  A((0  +  7r)g(oj  +  7t)  =  0, 


or  equivalently  in  matrix  notation. 


A(®)  ^(®  +  k) 

h((0)  g(co) 

— 

1  0 

^(®)  g(®  +  71 )_ 

h((0  +  Ti)  g(®  +  n)_ 

0  1. 

(2.6) 


Hence,  if  we  denote  modulation  matrices,  m((i))  and  m((o),  respectively  by 


w(co)  = 


/z(co)  A(co  +  n) 
Lg(co)  g(0)  +  7C)J 


and  m(&)  = 


h(Gi)  h((o  +  n) 
Lg((o)  g(0)  +  %)_ 


then,  the  biorthogonality  condition  (2.6)  becomes 

V®  e  R,  m(®)m'(ro)  =  I .  (2.7) 

By  transposing  equation  (2.6),  we  see  that 

;z(®)^(o)  +  g(®)g(®)  =1  ^2  8) 

/z(®)A(®  +  Tl) +  g(©)g(®  +  TT)  =  0. 

Equation  (2.8)  can  also  be  driven  by  a  general  two  channel  subband  coding  scheme  shown  in 
Figure  42.  A  discrete  input  signal  s(n)  is  decomposed  into  an  approximation  signal  a(n)  and  a 
detail  signal  d(n).  They  are  defined  by 
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Figure  42.  A  two  channel  subband  coding  scheme. 


fl(2co)  =  :^(A(C0)5(C0)  +  A((0  +  7t)5((0  +  7t))  , 

2  (2.9) 

dildi)  =  ^(g(C))5(CO)+g(G)  +  7t)5(Q)  +  7l))  . 

The  reconstructed  signal  in  the  frequency  domain  can  be  written  by 

r((o)  =  a(a))5'((o)  +  P(co)s'((0  +  7t)  ,  (2.10) 

where 

a(co)  =  /?(co)A(o))+^(co)g((o) , 

P(co)  =  /!(co)A(co  +  7r)  +  g(co)g(co  +  7t)  . 

In  order  to  achieve  perfect  reconstruction  in  the  subband  coding  system,  we  keep  only  the  linear 
shift-invariant  (LSI)  system  response  appearing  in  the  first  term  of  equation  (2.10)  and  cancel  out 
the  system  aliasing  in  the  second  term  of  equation  (2.10).  From  a(co)  =  1  and  P((o)  =  0,  we 
can  derive  equation  (2.8). 

By  applying  Cramer’s  rule  to  (2.8),  we  observe  that 

andg(<o)=-^^,  (2.12) 

DJe>)  0„(<») 

where  D^((o)  =  det  w((o) .  In  the  case  of  finite  filters,  we  need  to  restrict  D^(g>)  to  be  a 
monomial,  and  we  choose 

Z)Jco)  = (2.13) 
By  applying  equation  (2.13)  to  equations  (2.12)  and  (2.8),  we  have 
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g(a))  =  — e“'®/!((jo  +  7i)  and  g(co)  =  — +  ti)  , 


(2.14) 


and  equation  (2.8)  becomes 

h{G>)h{(o)  +  h{(ii  +  n)h{())  +  k)  =  1.  (2.15) 

The  filters  satisfying  equation  (2.15)  are  called  conjugate  filters. 

After  finding  the  conjugate  filters,  we  can  construct  a  scaling  function.  By  the  refinement 
relation  (2.1),  the  Fourier  transform  of  the  scaling  function  is  then  defined  by 

((.(to)  =  (21«) 

where  H  is  the  Fourier  transform  of  the  finite  filter  h,  in  other  words 

00 

/f((o)  =  ^  A(n)e“*”®. 

n  =  -00 


Since  (p(0)  =  1 ,  infinite  iterations  of  (2.16)  with  a  conjugate  filter // satisfying  (2.15)  yields  the 
limiting  scaling  function 

00 

9(0)  = 

7  =  1 

A  similar  statement  holds  for  the  dual  scaling  function  9  and  its  associated  dual  filter  H. 
Since  the  properties  of  the  filter  H  determine  the  properties  of  the  scaling  function  such  as  the 
smoothness  and  its  asymtotic  decay  at  infinity,  we  do  not  need  the  scaling  function  itself  in  many 
applications,  rather  we  directly  work  with  the  finite  filter  H. 

Figure  43  shows  the  limiting  process  to  make  the  Daubechies’s  four  tab  scaling  function. 
Each  figure  shown  in  Figure  43  shows  scaling  functions  without  iteration,  after  one  iteration,  after 
•  two  iterations  and  after  fifteen  iterations,  respectively,  from  top  to  bottom.  Figure  44  displays  a 
graph  of  mean  square  errors  involved  in  the  limiting  steps  to  make  the  Daubecies’s  four  tap 
scaling  function.  The  limiting  process  makes  little  difference  after  fifteen  iterations 

(RMS<2.24e-i6). 

As  pointed  out  in  [79],  useful  analyzing  properties  such  as  oscillations,  and  moments  can 
be  concentrated  on  the  vj;  whereas  the  synthesis  properties  such  as  regularity  are  determined  by 
\\f .  9  and  9  can  have  very  different  regularity  properties,  however  in  practice,  9  is  much  more 
regular  than  vj;  [16].  In  other  words,  analyzing  dual  wavelets  have  more  vanishing  moments  than 
synthesizing  wavelets. 
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Figure  43.  Graph  showing  the  limiting  process  to  make  the  Daubechies  tap  4 
scaling  function.  From  top  to  bottom,  scaling  functions  without  iteration,  after 
one  iteration,  after  two  iterations,  and  after  fifteen  iterations,  respectively. 


Figure  44.  Graph  of  root  mean  square  errors  involved  in  the  limiting  steps  to 
make  the  Daubechies  tap  4  scaling  function  shown  in  Figure  43. 
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2.4.4  Lifting  Scheme 


We  recall  the  Vetterli-Herley  lemma  [80]  which  defines  a  complete  characterization  of 
filters,  h  and  h  ,  biorthogonal  to  a  given  filter,  h.  The  two  dual  filters  h  and  h  are 
biorthogonal  to  h  in  the  sense  of  (2.15),  and  are  related  by 

h{Gi)  =  +  7r)5(2co)  ,  (2.18) 


where  5(00)  is  a  trigonometric  polynomial.  In  other  words,  if  one  of  the  dual  filters  is 
biorthogonal  to  h,  and  they  are  related  through  (2.18),  the  other  dual  is  biorthogonal  to  h  as  well 
[75].  By  combining  the  biorthogonal  condition  for  basis  functions  (2.4)  and  (2.18),  we  can  show 
the  following  relation  for  the  associated  FIR  filters.  From  an  initial  set  of  finite  biorthogonal  filters 
{ h,  g} ,  a  new  set  of  finite  biorthogonal  filters  can  be  found  by 


'V  old  . . — 

A(co)  =  h  {(i>)  +  g(G))s(2(o) 
g((o)  =  g®^'^(co)-^((o)j(2(o) . 


(2.19) 


This  is  called  lifting  and  this  operation  leads  to  the  lifting  scheme.  Primal  and  dual  basis  functions 
with  lifting  are  then  given  by: 

(p(x)  = 


(p(x)  =  ^(2x-k)  +  Y,S-k^(x-k) 

k  k 

\|/(x)  =  ^p''^‘^(x)-'^S|^(f>(x-k) 
k 

V(^)  =  25]gf‘^9(2^-^)  . 

k 


(2.20) 


There  are  several  facts  notable  in  (2.20).  The  scaling  function  does  not  change  after  lifting. 
The  wavelet  after  lifting  uses  the  scaling  function  at  the  same  level.  The  wavelet  construction 
formula  tells  us  that  the  property  of  the  wavelet  after  lifting  is  mostly  determined  by  the  property 
of  the  trigonometric  polynomial  s.  The  dual  wavelet  also  changes.  However  since  the  changes  in 
the  dual  wavelets  are  from  the  changes  in  the  dual  scaling  functions,  the  changes  of  the  dual 
wavelet  are  not  meaningful  by  themselves.  To  summarize,  the  lifting  algorithm  contributes  to 
some  fundamental  and  meaningful  changes  in  the  wavelet  and  the  dual  scaling  functions  through 
5  which  defines  the  true  power  of  the  lifting  scheme.  In  other  words,  a  multiresolution  analysis 
can  be  achieved  by  starting  from  a  simple  basis  function  and  building  a  more  performant  one  with 
desirable  properties. 

We  can  also  define  a  dual  lifting  from  an  initial  set  of  filters  { h,  g,  g°  }  and  a  dual 
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lifting  polynomial  s  such  that 

h((o)  =  h^^^((o)  +  g(oi))s(2co) 

g(G))  =  g'^^‘^((o)-h((o)s(2(o) . 

The  dual  lifting  scheme  is  thus  defined  by  simply  toggling  the  tildes  in  (2.20)  as 

^(x)  = 

(p(x)  =  2^11^"^  (p(2x-k)  +  J^S_^yif(x-k) 
k  k 

k 

¥(^)  =  22]gf'^(p(2x-A:)  . 

A: 


(2.21) 


(2.22) 


By  alternating  the  primal  and  dual  lifting  scheme,  one  can  come  up  with  a  multiresolution 
analysis  geared  for  specific  properties.  This  so-called  cakewalk  algorithm  leads  to  a  canonical 
form  of  the  wavelet  transform  shown  in  Figure  45.  Here  the  filters  followed  by  successive 
subsampling  define  the  Lazy  wavelet  transform  which  simply  splits  an  original  signal  into  odd  and 
even  subsets.  The  dual  lifting  predicts  wavelet  coefficients  with  the  help  of  the  scaling  function 
coefficients  based  upon  the  correlation  present  in  the  finer  signal.  The  primal  lifting  then  updates 
the  scaling  function  coefficients  using  the  previously  predicted  wavelet  coefficients  to  remove 

possible  aliasing  resulting  from  downsampling. 

Based  upon  the  formulae  of  the  primal  and  dual  lifting  schemes,  a  fast  forward  and  inverse 


Figure  45.  A  1-D  discrete  wavelet  transform  with  lifting  polynomials. 
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wavelet  transform  can  thus  be  formulated  as 


FORWARD  TRANSFORM 

(I)  Xji=  and  yj,i= 

k  k 

N-\ 

(II)  yj,i  =  jjj-'Zh-khk 

k 

N-\ 

(III)  Xjj  =  , 

k 


INVERSE  TRANSFORM 

^-1 

(I)  ijj  =  Xjj-  Yj  ^i-ky Ik 
k 

N-l 

m  yjj  =  yj.i*  S 

k 

(III)  X  .l  i  =  J2'£ht_2,Xjj+  . 

I  I 


Here  N  and  N  are  the  orders  of  the  primal  and  dual  lifting  polynomials  defined  as  s  and  5 
respectively. 

Now  we  can  define  a  canonical  implementation  structure  of  the  lifting  scheme  with  three 
stages:  split,  predict  and  update,  which  is  well  described  in  the  cakewalk  algorithm  shown  in 
Figure  45.  A  simpler  version  of  the  forward  wavelet  transform  with  lifting  is  shown  in  Figure  46. 
The  split  process  splits  the  input  data  { +  j }  into  two  smaller  subsets  { Xj }  and  { } .  From  the 
perfect  reconstruction  property  in  the  polyphase  representation,  the  most  trivial  subsets  are  even 
and  odd  samples.  This  is  often  referred  to  as  the  Lazy  wavelet  transform  [75]  as  described  earlier. 
In  the  second  stage  Predict,  we  try  to  use  the  { Xj }  subset  to  predict  the  { y^ }  subset  based  on  the 
correlation  present  in  the  original  data.  The  { Xj)  subset  is  then  lifted  with  the  help  of  the  wavelet 
coefficients  {y^} .  The  idea  is  to  find  a  better  }  so  that  a  certain  scalar  quantity  Q[),  such  as 
the  mean  is  preserved,  QfXj^f)  =  OSXj) .  These  aspects  lead  to  the  following  forward  wavelet 
transform: 
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yj 
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Figure  46.  Canonical  structure  in  implementing  1-D  discrete  forward  wavelet 
transform  with  lifting. 


Fory  =  -1  downto-n:  s 


=  sputCkj^^) 
Xj  =  Xj+uiyj) , 


where  fPand  U  are  prediction  and  update  operators,  respectively.  The  inverse  transform  is 
immediately  driven  by  reversing  the  operations  and  toggling  the  signs: 


For  j  =  —n  upto  —  1 :  "I 


yj 

^y+i 


y  j 

JoiniXj,  yj) 


Figure  47  shows  the  wavelet  representation  with  lifting  after  level  1  decomposition.  The 
wavelet  transform  presented  here  in  fact  is  the  (N  =  2,  N  =  2)  biorthogonal  wavelet  transform  of 
Cohen-Daubechies-Feauveau  [79]. 


2.4.5  Interpolating  Scaling  Functions  and  Associated  Wavelets 

Sweldens  [75]  [76]  [81]  showed  that  the  lifting  scheme  described  in  the  previous  section  is 
connected  with  interpolating  scaling  functions.  The  use  of  interpolating  scaling  functions  in 
multiresolution  analysis  has  been  studied  by  several  authors  [82]  [83].  The  main  advantage  of  the 
interpolating  scaling  function  is  that  the  coefficients  in  a  function  expansion 

fix)  =  '^X^ipix-k) 
k 

can  be  easily  obtained  from  =  /(k)  since  the  interpolating  scaling  function  is  defined  as 
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(a) 


(b) 


Figure  47.  Wavelet  representation  with  lifting  after  level  1  decomposition,  (a),  (b),  (c) 
and  (d)  show  DC  component,  wavelet  coefficients  processed  in  the  vertical  direction, 
in  the  horizontal  direction  and  in  the  diagonal  direction,  respectively. 
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(p(A:)  =  6^  for  all  A:  e  Z.  By  combining  the  definition  of  the  interpolating  scaling  fimction  and 
the  refinement  relation  (2.1),  we  get  h2f^  -  S^/2  which  becomes 

A(o))  +  A((o  +  7i)  =  1.  (2.23) 

Filters  satisfying  (2.23)  is  called  interpolating  filters,  a-trous  filters  are  the  interpolating  filters 
which  are  used  to  compute  coefficients  of  a  continuous  wavelet  transform  [84]  [66]. 

There  have  been  several  approaches  to  compute  interpolating  scaling  functions.  Here  we 
briefly  discuss  one  proposed  by  Deslauriers  and  Dubuc  [82]:  interpolating  subdivision,  which  is 
connected  to  our  Overcomplete  Lifting  Scheme.  The  interpolating  subdivision  is  formally  defined 
as  follows.  We  first  construct  a  polynomial  P  of  degree  N-l  so  that  Pixj^  j^  +  i)  =  ^j,k+i 
-D+l<l<D,  where  N  =  2D .  We  then  calculate  coefficients  on  the  next  finer  level  as  the 
value  of  this  polynomial  by  + 1,2^+ 1  =  P(^j+i,2kv  1)  •  ^8  shows  the  scaling  functions 

with  interpolating  subdivision  of  order  2, 4, 6,  and  8  from  top  to  bottom. 

These  interpolating  scaling  functions  have  compact  supports,  and  are  symmetric,  so  that  a 
(p(x)  is  zero  outside  the  interval  [-W+  1,  N- 1] .  The  subdivision  interpolating  scaling 
functions  and  their  translates  reproduce  polynomials  up  to  degree  -  1 .  They  also  satisfy  the 


-2  0  2 


-4  -2  0  2  4 

Figure  48.  Interpolating  scaling  functions  resulted  from  the  interpolating 
subdivision  of  order  2, 4, 6  and  8  fi*om  top  to  bottom. 
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refinement  equation  (2.1),  which  is  critical  in  defining  a  multiresolution  analysis.  Since  the  order 
of  the  reproduced  polynomials  is  invariant  within  an  interval,  the  resulting  scaling  functions  are 
well  adapted  on  an  interval  of  finite  signal  without  extensions  such  as  zero  padding,  periodization 
or  reflection.  Figure  49  shows  an  example  of  cubic  (N  =  4)  scaling  functions  changing  their  shapes 
to  accommodate  themselves  near  the  left  boundary  without  degrading  the  regularity.  The 
interpolation  properties  previously  described  are  still  preserved  well  on  the  boundary  for  this 
scaling  function. 

Table  2-1  lists  the  lifting  coefficients  needed  for  the  interpolating  subdivision  scheme 
which  are  computed  from  Neville's  algorithm.  Table  2- 1(a),  (b),  (c)  and  (d)  show  the  filter 
coefficients  ^oxN=2,N=A,N-6  and  JV  =  8,  respectively.  Each  lists  all  coefficients  to  compute 
values  at  every  possible  location  newly  introduced  as  a  middle  point  at  the  next  finer  level. 

By  filling  in  the  refinement  relation  in  the  interpolating  case,  we  see  that 


Figure  49.  The  cubic  {N  =  4)  interpolating  scaling  functions  at  k  =  1,  2,  3,  4,  and  5 
(top  to  bottom)  who  change  their  shapes  to  alleviate  boundary  effects. 
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#  of  X's  on 
(left,  right) 

Lifting  Coefficients 
(N  =  4) 

(0,4) 

(2.1875,  -2.1875,  1.3125,  -0.3125) 

(1,3) 

(0.3125, 0.9375,  -0.3125,  0.0625) 

(2,2) 

(-0.0625, 0.5625,  0.5625,  -0.0625) 

(3,  1) 

(0.0625,  -0.3125, 0.9375, 0.3125) 

(4,0) 

(-0.3125, 1.3125,  -2.1875, 2.1875) 

(b) 


UofX's  on 
(left,  right) 

Lifting  Coefficients 
(N  =  2) 

(0,2) 

(-0.5,  1,5) 

(1,1) 

(0.5,  0.5) 

(2,0) 

(1.5, -0.5) 

#  of  X's  on 
(left,  right) 

Lifting  Coefficients 
(N  =  6) 

(0,6) 

(2.7070,  -4.5117,  5.4141,  -3.8672, 1.5039,  -0.2461) 

(1,5) 

(0.2461, 1.2305,  -0.8203, 0.4922,  -0.1758, 0.0273) 

(2,4) 

(-0.0273, 0.4102, 0.8203,  -0.2734, 0.0820,  -0.01 17) 

(3,3) 

(0.0117,  -0.0977, 0.5859, 0.5859,  -0.0977, 0.0117) 

(4,2) 

(-0.01 17, 0.0820,  -0.2734, 0.8203, 0.4102,  -0.0273) 

(5,1) 

(0.0273,  -0.1758, 0.4922,  -0.8203, 1.2305, 0.2461) 

(6,0) 

(-0.2461, 1.5039,  -3.8672,  5.4141,  -4.5117,  2.7070) 

(c) 


#  of  X's  on 
(left,  right) 

Lifting  Coefficients 
(N  =  8) 

(0,  8) 

(3.1421,  -7.3315, 13.1968,  -15.7104, 12.2192,  -5.9985, 1.6919,  -0.2095) 

(1,7) 

(0.2095, 1.4663,  -1.4663, 1.4663,  -1.0474, 0.4888,  -0.1333, 0.0161) 

(2,6) 

(-0.0161, 0.3384, 1.0151,  -0.5640, 0.3384,  -0.1450, 0.0376,  -0.0044) 

(3,  5) 

(0.0044,  -0.0513, 0.4614, 0.7690,  -0.2563, 0.0923,  -0.0220, 0.0024) 

(4,4) 

(-0.0024, 0.0239,  -0.1196,0.5981, 0.5981,  -0.1196, 0.0239,  -0.0024) 

(5,  3) 

(0.0024,  -0.0220, 0.0923,  -0.2563, 0.7690, 0.4614,  -0.0513,  0.0044) 

(6,2) 

(-0.0044,  0.0376,  -0.1450, 0.3384,  -0.5640, 1.0151, 0.3384,  -0.0161) 

(7,  1) 

(0.0161,  -0.1333, 0.4888,  -1.0474, 1.4663,  -1.4663,  1.4663, 0.2095) 

(8,0) 

(-0.2095, 1.6919,  -5.9985, 12.2192,-15.7104,  13.1968,  -7.3315,3.1421) 

(d) 


Table  2-1.  Lifting  coefficients  based  on  the  number  of  X's  on  its  left  and  right,  (a),  (b), 
(c)  and  (d)  show  the  filter  coefficients  for  7^=  2,  iV=  4,  A^=  6  and  iV=  8,  respectively. 
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^k,K  ~  ~  ^^j,k,l^j+\,l^^j+\,2k)  ~  ^j,k,2k' ■> 

which  can  be  written  as 

''^j,k  ~  ^j+\,2k'^'^^j,k,2m  +  \'^j+\,2m  +  \  • 
m 

An  approximation  Ajf  from  a  finer  one  Aj+  j/  is  then  built  by  simple  subsampling, 

Xj  =  Xj+i^2k’  followed  by  one  step  subdivision  on  Aj : 

jf  ~  “  ^^j+l,2k^j+l,2k'^^^^j,k^j,k,2m  +  l^J+^,2m  +  l 

^  K  Km 

The  detail,  + 1  —  ,  consists  only  of  fimctions  of  the  form  tPy  + 1, 2m  + 1  •  ^1^^^ 

wavelets  are  given  by  \|/y  ^  =  (Py  + 1  2m  + 1  •  However,  a  necessary  condition  for  the  wavelets  to 
form  a  Riesz  basis  is  that  the  dual  vanishing  moment  N  is  at  least  one.  A  new  wavelet  with  N  is 
then  formulated  by  following  with  lifting  scheme: 

N-\  j 

^J,m  ~  ^y+ l,2m  + 1  ~  m^y, m  +  /  ’ 

1  =  0 

where  the  lifting  coefficients  Sj^  ^  are  chosen  so  that 

f'Ky.mW  =  0,  JxVj/y  ,„(x)  =  0,  . .  Jx^->V/y,mW  =  0  . 

In  other  words,  the  new  wavelet  is  composed  of  the  old  wavelet,  a  subsampled  scaling 
function  located  at  odd  indices,  and  its  immediate  neighboring  scaling  functions  at  the  same  scale 
j.  Figure  50  shows  wavelets  with  two  dual  vanishing  moments,  i.e.  N  =  2,  associated  with 
interpolating  scaling  functions  the  order  of  the  subdivision  of  which  are  2, 4, 6,  and  8  from  top  to 
bottom.  Figure  51  shows  a  wavelet  accommodating  the  left  boundary.  The  cubic  (N=  4)  scaling 
function  was  used  to  construct  this  wavelet  with  N  =  2  vanishing  moments. 

2.4.6  Generalized  Lifting  Scheme 

We  have  seen  in  the  previous  sections  that  there  are  many  advantages  of  the  lifting  such  as 
an  in-place  implementation  of  the  fast  wavelet  transform,  no  need  for  the  Fourier  transform, 
easiness  of  the  inverse  transform,  capability  of  custom-design  of  wavelets  and  adaptiveness  of 
wavelet  transforms.  However  the  lifting  scheme  discussed  was  used  in  constructing  biorthogonal 
wavelets.  However  can  the  lifting  be  used  in  constructing  any  wavelets  or  wavelet  transforms  with 
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Figure  50.  Wavelets  with  two  dual  vanishing  moments  {N  =  2)  associated  with 
interpolating  scaling  functions  whose  orders  of  the  subdivision  are  2,  4,  6,  and  8  from 
top  to  bottom. 

finite  filters?  Daubechies  and  Sweldens  addressed  the  generality  of  the  lifting  by  connecting  the 
Euclidean  algorithm  to  a  subband  transform  in  the  z-domain  [77]. 

The  z-transform  of  a  FIR  filter  his  a  Laurent  polynomial  H(z)  given  by 

b 

=  Z 

k-  a 

where  a  and  b  are  respectively  the  lower  and  upper  bound  of  the  support  of  the  filter  h.  The  degree 
of  the  Laurent  polynomial  H{z)  is  given  by  \H{z)\  =  a,  so  that  a  polynomial  Cz  has  zero 
degree  and  is  a  monomial,  where  C  is  a  non-zero  constant  and  /  is  an  integer.  A  Laurent 
polynomial  is  invertible  if  and  only  if  it  is  a  monomial. 

The  monomial  constraint  plays  a  key  role  in  finding  a  modulation  matrix  or  a  polyphase 
matrix  in  a  filter  bank  coding  system.  Figure  52  shows  a  discrete  wavelet  transform  in  z-domain. 
The  conditions  for  perfect  reconstruction  are  given  by 


113 


i/(z)^(z-i)  +  G(z)G(z-i)  =  2  , 
i/(z)^(-z-i)+  c;(z)d(-z-i)  =  0  . 


(2.24) 


The  modulation  matrices  M(z)  and  M(z)  are  defined  as 

M(z)  =  ,  M(z)  =  .  (2.25) 

_G(z)  G(-z)J  LG(z)  G(-z)_ 

By  using  these  modulation  matrices,  the  perfect  reconstruction  condition  (2.24)  can  be  written  as 

M(-z-1)^M(z)  =  21,  (2.26) 

where  I  is  the  2  by  2  identity  matrix. 

We  can  simplify  the  perfect  reconstruction  in  the  filter  bank  structure  using  a  polyphase 
representation  which  provides  an  efficient  and  convenient  tool.  The  polyphase  implementation 
works  on  different  phases  separately.  The  input  is  separated  into  even  and  odd  phases  followed  by 
subsampling.  Then  the  separate  phases  of  the  filters  act  in  parallel  on  the  separate  phases  of  the 
input.  Figure  53  shows  a  polyphase  representation  of  the  wavelet  transform  shown  in  Figure  52. 
The  polyphase  representation  of  a  filter  H  is  given  as  H{z)  =  Hg(z^)  +  z~^H^(z^) ,  where 
and  contain  respectively  the  even  and  odd  coefficients.  If  we  denote  the  polyphase  matrix 
P(z)  as 

(2.27) 

then  the  relation  between  the  modulation  matrix  (2.25)  and  the  polyphase  matrix  (2.27)  becomes 


LP 

HP 

Figure  53.  Polyphase  representation  of  the  wavelet  transform. 


/>(.)  =  . 
H,iz)  C?„(z) 
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(2.28) 


P(z^) 


\  ^  ^  MKz) 


1 

2 


M{z) 


Now  the  perfect  reconstruction  property  with  modulation  matrices  (2.26)  can  be  written  in 
polyphase  matrices  as 

P{z)P{z-^y  =  I .  (2.29) 


From  the  invertible  property  of  Laurent  polynomials  and  (2.29),  the  determinant  of  P{z)  must  be 
a  monomial  in  z  so  that  P(z)  =  Cz^ .  In  other  words,  finding  FIR  filters  for  a  wavelet  transform 
becomes  finding  a  polyphase  matrix  P(z)  with  the  determinant  one.  By  combining  this  monomial 
property  and  (2.29),  the  dual  polyphase  matrix  P{z)  is  given  as 


P(z)  = 


The  most  trivial  case  of  the  polyphase  matrix  satisfying  the  condition  of  perfect 
reconstruction  is  P(z)  =  I  which  implies  ff(z)  =  H{z)  =  land  G(z)  =  G(z)  =  z“* .  The 
transform  with  P(z)  =  I  is  often  referred  as  the  Lazy  wavelet  transform  since  it  does  nothing  but 
separating  the  input  into  even  and  odd  subsets. 

Using  the  polyphase  structure,  the  lifting  given  in  (2.19)  at  page  104  can  be  written  in  z- 
domain  as 

H(z)  =  +  ^)5(z-2)  ^2.30) 

G(z)  =  G^‘‘^(z)-Hiz)S(z^) , 

and  the  dual  lifting  in  (2.21)  at  page  105  as 

H(z)  =  Ho^df^z)  +  ^)5(z-2) 

G(z)  =  G^^‘^(z)-H(z)~S(z^)  . 

There  are  two  noble  identities  for  a  decimator  and  interpolator.  In  most  applications  a 
decimator  is  preceded  by  a  bandlimiting  filter  H{z)  to  avoid  aliasing.  An  interpolation  filter,  on 
the  other  hand,  is  followed  by  an  interpolator  to  eliminate  imaging  effect.  Figure  54  shows 
equivalent  identity  systems.  In  Figure  54(a),  we  have  a  decimator  followed  by  a  transfer  function 
F(z) .  This  cascade  is  equivalent  to  the  one  in  the  right  provided  F(z)  is  a  rational  transfer 
function,  i.e.  a  ratio  of  polynomials  in  .  In  Figure  54(b),  we  show  that  the  two  cascades  in  the 
left  and  right  are  equivalent  for  a  interpolator.  Based  on  the  property  of  the  noble  identities  with 
M  =  2  and  L  =  2 ,  the  polyphase  components  of  H(z)S(z^)  are  Hg(z)S{z)  for  even  and 
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Figure  54.  Noble  identities  for  multirate  systems.  Identity  systems  for  (a)  a  decimator 
and  (b)  an  interpolator. 
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where  is  the  GCD,  n  is  the  smallest  number  for  which  B^{z)  =  0,  and 

Qf{z)  =  A^_^{z)/B|_^{z) .  Here  if  A^(z)  is  a  monomial,  then  A(z)  and  B{z)  are  relatively 
prime.  We  can  always  choose  the  quotients  Qi(z)  so  that  the  GCD  A^(z)  is  a  constant. 

By  applying  the  Euclidean  algorithm  to  polyphase  matrices,  we  can  factor  any  pair  of 
complementary  filters  (H(z),  G(z))  into  lifting  steps.  A  filter  pair  (H(z),  G(z))  is  defined  to  be 
complementary  if  the  determinant  of  the  corresponding  polyphase  matrix  P(z)  is  one.  The 
polyphase  matrix  can  be  given  as 


P(z)  = 


m 

n 

i=  1 


-5,(z) 

1 


-lA 


1  0 
Siiz)  1 


K  0 

ok 


(2.37) 


and  the  dual  polyphase  matrix  as 
P(z) 


m 

n 

V  =  1 


1  0 
5,(z-i)  1 


1 

0  1 


1 


0 


K 
0  K 


(2.38) 


Tt 

Here  AT  is  a  non-zero  constant,  and  w  =  +  1  •  The  polyphase  matrix  for  analysis  filter  bank 

^  jL 

P{z~^y  is  thus  written  as 


P(z-i) 


-G?  = 


ko 

0  K 


1 

n 

M  =  W 


1  0 
-Siiz)  1 


1  ^.(z) 

0  1 


(2.39) 


Figure  55  shows  the  forward  wavelet  transform  using  the  dual  polyphase  matrix  with  finite 


Figure  55.  Forward  wavelet  transform  using  a  dual  polyphase  matrix  P(z“*)^with 
finite  number  of  lifting  steps. 
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Figure  56.  Inverse  wavelet  transform  using  the  polyphase  matrix  P{z)  with  finite 
number  of  lifting  steps. 


lifting  steps  given  in  (2.39).  Figure  56  represents  the  finite  steps  of  the  inverse  wavelet  transform 
using  the  polyphase  matrix  representation  in  (2.37). 

Every  wavelet  or  wavelet  transform  can  thus  be  obtained  with  a  factorization  of  the 
polyphase  matrix  representation  of  the  filter  bank  by  means  of  Euclidean  algorithm. 

2.4.7  Overcomplete  lifting  scheme 

In  this  section  we  propose  an  overcomplete  lifting  scheme  which  provides  a  framework  for 
multiscale  analysis  with  overcomplete  wavelet  representations.  As  described  in  the  previous 
chapter,  the  lifting  scheme  is  a  flexible  tool  for  constructing  compactly  supported  second 
generation  wavelets  which  are  not  necessarily  translates  and  dilates  of  one  mother  wavelet  [75]. 
We  have  seen  that  the  lifting  scheme  uses  a  simple  basis  function  and  builds  a  more  performant 
one  with  desirable  properties.  Flexibility  afforded  by  the  lifting  scheme  allows  basis  functions  to 
change  their  shapes  near  the  boundaries  without  degrading  regularities.  The  lifting  scheme  also 
provides  faster  processing  by  making  optimal  use  of  similarities  between  high  and  low  pass 
filters.  The  lifting  algorithm  was  originally  introduced  to  construct  biorthogonal  wavelets 
associated  with  interpolating  scaling  functions.  Daubechies  and  Sweldens  later  generalized  it  to 
construct  any  wavelets  with  a  finite  number  of  lifting  steps  [77]. 

The  lifting  scheme  utilizes  a  elassical  two-channel  filter  bank  as  a  framework  for 
multiresolution  analysis.  However  this  traditional  framework  is  not  translation  invariant  [15]. 
Representations  with  a  translation-invariant  characteristie  are  highly  desirable  for  feature 
analysis.  Mallat  and  Zhong  introdueed  an  adaptive  sampling  based  upon  local  maxima  as  an 
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overcomplete  wavelet  representation  [15]. 

In  this  section  we  address  the  following  question:  Can  the  lifting  scheme  be  used  as  a 
framework  for  overcomplete  wavelet  representations  with  feature  analysis  in  mind?  Our  work 
addresses  this  question  by  investigating  each  stage  of  the  multiscale  analysis  described  in  [75], 
We  use  a  smoothing  Lazy  wavelet  in  the  split  stage  which  does  not  subsample,  but  smooths  an 
input  image  so  that  the  low-pass  channel  contains  some  redundant  information.  We  then  show  that 
the  following  predict  and  update  stages,  based  upon  specific  properties,  indeed  lead  to  a  useful 
multiresolution  analysis. 

We  showed  how  the  (primal)  lifting  equation  (2.20)  and  the  dual  lifting  equation  (2.22) 
were  driven  and  how  they  can  be  embedded  in  a  framework  for  a  multiscale  analysis  which  is 
implemented  in  the  filter  bank  coding  system  shown  in  Figure  45.  This  computational  structure 
was  simplified  later  in  a  semantic  form  in  Figure  46.  Here  the  filters  followed  by  successive 
subsampling  was  defined  as  the  Lazy  wavelet  transform  which  simply  splits  an  original  signal  into 
odd  and  even  subsets.  The  dual  lifting  predicted  wavelet  coefficients  with  the  help  of  the  scaling 
function  coefficients  based  upon  the  correlation  present  in  the  finer  signal.  The  primal  lifting  then 
updated  the  scaling  function  coefficients  using  the  previously  predicted  wavelet  coefficients  to 
remove  possible  aliasing  resulted  from  the  downsampling. 

For  overcomplete  representations,  we  introduce  a  smoothing  lazy  wavelet  in  the  split  stage 
which  does  not  downsample  the  input  image  so  that  + 1,  *  ="  “  Yy,  m  ■  Here  we  choose  the 

indices  so  that  j  eJ,  I  e  K(j) ,  k  €  K(j  +  1) ,  and  m  e  M{j)  where  M(J)  =  K(j  +  1  )\  K(j) . 
We  then  smooth  ^  using  a  smoothing  filter,  which  is  often  a  Gaussian  [15].  This  new  splitting 
methodology  in  the  filter  bank  scheme  is  then  followed  by  the  decorrelation  of  wavelet 
coefficients  using  the  dual  lifting.  Since  the  new  split  operation  does  not  use  downsampling  to 
make  distinct  subsets  { ,}  and  {jj^ ,} ,  there  is  no  aliasing.  The  primal  lifting  introduced  to 
compensate  for  the  aliasing  by  preserving  energy  between  two  contiguous  approximations  needs 
no  longer  exist.  The  new  forward  and  inverse  wavelet  transform  algorithm  for  overcomplete 
representations  can  thus  be  formulated  as 


FORWARD  TRANSFORM 

bplit.  Xj  i  —  hk—2iX j  + 1  1  j,m  ~  ^y+i,A 

k 

Dual  Lifting:  y j  , 


(2.40) 


.m  —  ln 


n  =  0 
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INVERSE  TRANSFORM 


N-l 

[nverse  Dual  Lifting:  y j  „,  = 

n  =  0 

^oin:  =  Xj  i  +  yj  ,„ 

where  K(j  +  1 ) ,  L(j)  and  M{j)  are  index  sets  such  that  \K(J  +  1)|  =  \M (y)!  =  l■^(7)l » 

{k\k  e  K(j  +  1)} ,  {l\leL(j)}  andm\m  g  M(j).  N  and  are  primal  and  dual  vanishing 
moments,  respectively.  Figure  57  shows  an  undecimated  version  of  two  channel  analysis/ 
synthesis  filter  bank  with  the  newly  proposed  overcomplete  lifting  scheme.  By  iterating  this 
scheme  through  the  low-pass  channel,  we  can  achieve  the  1-D  discrete  wavelet  transform  for 
overcomplete  multiscale  representations. 

In  the  previous  section  we  discussed  how  interpolating  scaling  functions  are  related  to  the 
lifting  scheme  through  the  interpolating  subdivision  proposed  by  Deslauriers  and  Dubuc  [82]. 
Here  we  use  the  same  mechanism,  but  the  approach  is  different.  Instead  of  refining  samples  by 
filling  each  middle  point  between  two  contiguous  samples  at  a  finer  level  using  subdivision  with  a 
constant  spatial  support,  we  keep  the  size  of  the  samples  constant  throughout  the  decomposition, 
but  change  the  support  of  the  filter  associated  with  the  interpolating  scaling  function. 

Figure  58  illustrates  how  to  compute  the  wavelet  coefficients  with  a  cubic  dual  lifting 
polynomial  {N=  4).  In  this  figure,  we  show  how  the  samples  at  the  next  analyzing  level  are  related 
to  the  samples  at  the  finer  level  using  two  points  one  of  which  is  near  the  left  boundary,  and  other 


■In 


(2.41) 


Figure  57.  Undecimated  version  of  a  two  channel  filter  bank  with  the  proposed 
overcomplete  lifting. 
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Figure  58.  Forward  wavelet  transform  to  extract  detail  information  as  failures  to  be 
cubic. 


not.  The  predicted  wavelet  coefficients  contain  the  amount  by  which  the  coarser  approximation 
locally  fails  to  be  cubic.  Here  the  dual  lifting  coefficients  s  =  (-0.0625, 0.5625, 0.5625,  -0.0625}. 
The  coefficients  are  adjusted  near  the  boundaries  of  the  finite  input  signal  to  avoid  boundary 
effects.  For  A^=  6,  5  =  (0.0117,  -0.0977, 0.5859, 0.5859,  -0.0977, 0.0117}  and  5  =  (-0.0024, 
0.0239,  -0.1196,  0.5981,  0.5981,  -0.1196,  0.0239,  -0.0024}  for  A^=  8. 

For  images,  2D  masks  can  be  used.  Figure  59(a)  is  a  mask  for  A/^=  4  which  is  used  for 
regions  which  are  not  affected  by  the  boxmdaries.  Figure  59(b)  is  for  the  pixel  (0,  0)  of  an  image. 
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Figure  59.  Examples  of  2-D  masks  of  dual  lifting  coefficients  for  N=  4.  (a)  A  mask  for 
regions  which  are  not  affected  by  the  boundaries,  (b)  A  mask  for  the  point  (0,  0)  of  an 
input  image.  Rectangles  indicate  points  that  are  modified  by  masking  operations. 
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The  proposed  algorithm  has  been  applied  to  a  mammogram  of  size  512  by  512.  We  used  a 
cubic  spline-based  filter  for  // described  in  [15]  for  the  smoothing  lazy  wavelet  transform.For  the 
dual  lifting,  we  applied  a  cubic  polynomial  with  coefficients  5  =  {-0.0625 , 0.5625 ,  0.5625 , 
-0.0625 }.  Figure  60  shows  wavelet  representations  for  three  levels  of  analysis.  We  used  the 
wavelet  coefficients  to  reconstruct  an  image  with  an  error  on  the  order  of  1 1.3dB.  Since  the  lifting 
scheme  can  speed  up  wavelet  processing  by  a  factor  of  two  [76]  and  convolution  using  a  block 
algorithm  can  ftxrther  reduce  processing  time  by  two-fold  within  each  analyzing  level  [85],  the 
proposed  algorithm  can  achieve  an  overall  speedup  of  4  over  the  standard  overcomplete  wavelet 
transform  where  an  FFT  is  the  basic  implementation  tool.  As  pointed  out  in  [85],  the  FFT  method 
performs  best  when  the  filter  length  is  long  (over  16  taps).  These  results  suggest  that  our  approach 
for  computing  overcomplete  wavelet  representations  with  compactly  supported  lifting  filters, 
compared  to  existing  algorithms,  is  significantly  more  efficient. 


2.4.8  Summary 

We  have  shown  that  the  lifting  algorithm  works  well  in  more  general  domain  since 
wavelets  computed  from  lifting  change  shape  near  boundaries  without  degrading  regularity.  This 
property  alleviates  boundary  artifacts  which  appear  in  existing  algorithms.  Thus  our  algorithm 
provides  a  robust  way  of  computing  wavelet  coefficients  of  mammograms  for  arbitrary  domains 

(ROI). 

We  have  shown  an  approach  for  computing  overcomplete  wavelet  representations  wi 
compactly  supported  lifting  filters.  The  Overcomplete  Lifting  Scheme  performed  more  efficiently 
than  traditional  (existing)  methods.  This  algorithm  will  provide  the  computational  speed  up, 
needed  to  carry  out  methods  of  multiscale  enhancement  processing  interactively,  on  diagnostic 
electronic  display  systems  in  the  future. 
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Figure  60.  Wavelet  coefficients  of  a  mammogram  with  the  proposed  scheme.  The  first  and 
second  column  show  the  wavelet  coefficients  in  horizontal  direction  and  vertical  direction, 
respectively. 


Original  Mammogram 

Reconstructed  Mammogram 
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3  Conclusions 

In  the  paragraphs  below,  we  summarize  the  results  and  progress  made  during  the  past  year 
and  identify  the  final  tasks  of  our  investigation  to  be  carried  out  during  the  last  year  of  the 
project. 

We  extended  the  one-dimensional  discrete  dyadic  wavelet  transform  to  two  dimensions  and 
showed  that  in  case  of  wavelets  being  equal  to  the  second  derivative  of  a  central  B-spline 
function,  the  two-dimensional  discrete  dyadic  wavelet  transform  decomposition  is  included 
in  the  multiscale  spline  derivative-based  transform  decomposition.  This  attests  to  the 
flexibility  of  the  latter  transform,  which  can  incorporate  a  variety  of  wavelet  based 
mammographic  image  enhancement  methods. 

The  multiscale  spline  derivative-based  transform  was  used  for  a  novel  image  enhancement 
scheme:  the  input  image  was  first  processed  for  enhancement  of  specific  mammographic 
features,  and  the  resultant  enhanced  images  were  fused  into  a  final  result.  The  described 
scheme  represents  a  synthesis  between  our  early  global  methods  developed  for  an  overall 
improvement  of  image  quality  and  the  local  processing,  targeting  specific  mammographic 
features. 

The  employed  transform  is  highly  redundant,  and  efficient  implementation  is  of  particular 
importance.  We  presented  a  fast  filter  bank  implementation  that  takes  advantage  of  filter 
and  signal  symmetries.  Tests  of  our  method  included  convincing  quantitative  evidence  in 
terms  of  improved  visibility  and  detection  of  subtle  features. 

We  derived  and  formulated  a  mass  detector,  based  on  multiscale  analysis,  that  very  finely 
samples  the  scale  parameter,  in  search  for  subtle  masses  in  mammography.  The  importance 
of  arbitrary  scale  analysis  was  demonstrated  by  digital  radiographs  of  a  mammography 
phantom  and  digitized  mammograms.  For  the  first  time,  wavelets  and  CFAR  detectors 
were  combined  into  a  single  formulation.  We  demonstrated  that  this  powerful  detector  was 
able  to  detect  very  subtle  masses,  which  were  rated  to  be  almost  invisible  by  radiologist 
specializing  in  mammography. 

We  embedded  a  maxima-across-scale  scheme  [51,  52]  under  a  wavelet  framework.  And 
showed  that  wavelet  maxima  across  scales  can  provide  boundary  information  of  an  object 
and  indicate  its  location.  By  combining  wavelet  theory  and  the  geometric  knowledge  of  a 
tumor,  we  propose  to  use  a  normalized  contrast  as  an  index  to  find  the  best  scale.  This 
method  appeared  to  be  promising  in  selecting  the  best  scale  (automatically)  for  tumor 
detection. 

An  enhancement  algorithm  relying  on  multiscale  wavelet  analysis  and  extracting  oriented 
information  at  each  scale  of  analysis  was  investigated.  The  evolution  of  wavelet  coefficients 
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across  scales  characterised  the  local  shape  of  irregular  structures.  Using  oriented 
information  to  detect  the  features  proved  to  be  an  effective  method  of  enhancing  complex 
and  subtle  structures  of  the  breast.  Steerable  filters,  rotated  at  arbitrary  orientations 
reliably  found  visual  cues  within  each  spatial-frequency  sub-band  of  an  image.  Coherence 
measure  and  dominant  orientation  clearly  helped  to  discriminate  features  from  complex 
surrounding  tissue  in  mammograms. 

The  lifting  scheme  alone  can  speed  up  wavelet  analysis  processing  by  a  factor  o  two 
In  addition,  since  convolution  using  a  block  algorithm  can  further  reduce  processing  time 
by  two  told  within  each  level  of  analysis  (85],  our  proposed  algorithm  can  achieve  an  overaU 
speedup  of  four  over  standard  overcomplete  wavelet  transform  computations  where  an  FF 
is  the  basic  implementation  tool.  As  pointed  out  in  (321,  the  FFT  method  performs  best 
when  the  filter  length  is  long  (over  16  taps).  These  results  suggest  that  our  approach  or 
computing  overcomplete  wavelet  representations  with  compactly  supported  lifting  filters, 
compared  to  existing  algorithms,  is  significantly  more  efficient  than  previously  known 
methods.  This  algorithm  will  provide  the  computational  speed  up,  needed  to  carry  out 
methods  of  multiscale  enhancement  processing  interactively,  on  diagnostic  electronic 

display  systems  in  the  future. 

As  mentioned  in  the  overview  of  this  report,  the  final  report  of  the  project  will  contain  e 
results  of  an  ongoing  ROC  study.  This  study  is  being  carried  out  in  collaboration  with 
radiologists  and  medical  physicists  at  Columbia-Presbyterian  Hospital  at  Columbia 

University. 
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